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ABSTRACT 


This thesis explores using a direct pseudospectral method for the solution of 
optimal control problems with mixed dynamics. An easy to use MATLAB optimization 
package known as DIDO is used to obtain the solutions. The modeling of both low thrust 
interplanetary trajectories as well as aerocapture trajectories is detailed and the solutions 
for low thrust minimum time and minimum fuel trajectories are explored with particular 
emphasis on verification of the optimality of the obtained solution. Optimal aerocpature 
trajectories are solved for rotating atmospheres over a range of arrival V-infinities. 
Solutions are obtained using various performance indexes including minimum fuel, 
minimum heat load, and minimum total aerocapture mass. Finally, the problem 
formulation and solutions for the mixed dynamic problem of low thrust trajectories with a 
terminal aerocapture maneuver is addressed yielding new trajectories maximizing the 
total scientific mass at arrival. 
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I. INTRODUCTION 


In recent years, two developing concepts have been considered by mission 
designers to reduce the required propellant mass for interplanetary missions, thus 
allowing for increased mass of scientific payloads. Low thrust propulsion allows for 
greatly improved fuel efficiency and has the potential to increase payload mass fractions 
as well as providing trajectories not possible with impulsive thrust. Aerocapture is the 
careful management of a hypersonic atmospheric pass prior to orbit insertion to remove 
excess arrival velocity of a spacecraft obviating the need for a large propulsive capture 
maneuver. While aerocapture greatly reduces the propellant mass required for orbit 
insertion, the large heat loads generated in the atmospheric portion of the trajectory 
require the addition of a possibly massive thermal protection system (TPS). Both low 
thrust and aerocapture trajectories are highly non-linear with no closed-form solutions for 
their optimal open-loop control histories. This thesis addresses the feasibility of 
simultaneously optimizing a low thrust interplanetary trajectory with a terminal 
aerocapture maneuver to maximize the final scientific payload mass. In addition, the 
characteristics of optimal low thrust and aerocapture trajectories are explored with 
particular emphasis on techniques for the verification of the optimality of the obtained 
solutions. 
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H. LOW THRUST MODELING 


A. COORDINATE SYSTEM 

Polar coordinates were chosen for the low thrust portion of the trajectory design 
optimization. With the exception of the planets Mercury and Pluto, the orbits of the 
remaining planets all lie within 3.4° inclination of the solar ecliptic. For the purposes of 
this thesis, the small angle approximation has been used to assume that all planets’ orbits 
lie in the solar ecliptic. Figure 1 depicts the coordinate system used for the low thrust 
portion. 
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Figure 1: Low Thrust Coordinate System 
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B. EQUATIONS OF MOTION 


The state vector of a vehicle conducting a generic thrusting orbit transfer can be 
given as x = [r ,0, v r , v t , m ] 1 where, 

r = radial position from central body 
0 = transfer angle 

v, = radial velocity component 

v t = transverse velocity component 

m = vehicle mass 

The equations of motion for a constant thrust orbit transfer are given in [Ref. 1] 
and can be modified for non-continuous thrust as follows: 



( 1 ) 


J0_ = v L 
dt r 


( 2 ) 


dv r _ v, 2 (i T sinri 
dt r r 1 m 


( 3 ) 


dv t _ v r v t T cost) 
dt r m 


( 4 ) 


dm _ T 
dt v 

where 


( 5 ) 


(l = gravitational constant of central body 

T = thrust magnitude 

T] = thrust direction angle 

v e = exhaust velocity 
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c. 


GRAVITATIONAL MODEL 


Only the gravitational effects of the central body are considered in this work. 
Higher fidelity models include the gravitational effects of other celestial bodies thus 
enabling the discovery of gravity assist trajectories. For the case of general planet-to- 
planet orbit transfer, the vehicle is far from other gravitational bodies for the bulk of the 
trajectory. By assuming that the trajectory begins and ends at the boundary of the sphere 
of influence (SOI) with an escape velocity greater than or equal to zero, the gravitational 
effects of the origin and target planet can be neglected completely. 


For the inner, less massive planets, the SOI are reasonably small when compared 
to the orbit radii. Eqn. (6) [Ref. 2] approximates SOI radius as a function of the mean 
orbit radius and the ratio of masses between the gravitating bodies. 


m 




planet 


m.. 


' sun—planet 


( 6 ) 


Note that r S0I approximately scales linearly with the distance to the sun. Despite the fact 

that the SOI radii of the other planets increase with distance from the sun, this increase is 
roughly proportional to the increase in trajectory length. Considering this, it is 
reasonable to assume that the origin and target planet can be thought of as non-gravitating 
point masses, reducing the problem to a two-body problem. Table 1 gives the 
approximate SOI radius for selected planets. 
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Planet 

SOI Radius (km) 

SOI Radius (AU) 

Mercury 

1.13x10 s 

7.55xl0~ 4 

Venus 

6.17x10 s 

4.12x10 3 

Earth 

9.24x10 s 

6.18x10 3 

Mars 

5.74x10 s 

3.84xl0“ 3 

Jupiter 

4.83xl0 7 

0.32 

Neptune 

8.67X10 7 

0.58 


Table 1: SOI Radius for Selected Planets (adapted from [Ref. 2]) 


D. LAUNCH VEHICLE MODEL 

A launch vehicle’s performance is generally measured by the amount of mass it 
can place in a given orbit around the origin planet. For interplanetary trajectories, we are 
concerned with the trade-off between launch mass and the square of the escape velocity 
(know cryptically as C 3 ). [Ref. 2]. The launch performance may be fitted to a curve to be 
used to estimate launch mass from Q and vice versa. For this work, all trajectories will 
begin at the zero C 3 point (i.e. maximum launch mass with escape energy, equivalent to a 
parabolic escape). 


E. PROPULSION MODEL 

The low thrust propulsion system used for this research is modeled after the 
NASA Solar electric propulsion Technology Applications Readiness project (NSTAR). 
This engine, employed in the Deep Space I (DS1) mission [Ref. 4], uses a single xenon 
ion propulsion system. This engine is 30 cm in diameter with a mass of approximately 8 
kg. The engine delivers approximately 92 mN of thrust with a specific impulse of 3300 
seconds. To achieve higher thrust levels, multiple engines can be used. The engine’s 
specific impulse varies with the output thrust which is itself a function of the electrical 

power delivered to the engine. This model was simplified to assume a constant specific 
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impulse over the entire thrust range. Moreover, no limitations were placed on the 
available power provided to the engine. This assumption is valid provided a “near¬ 
limitless” source of power is available such as Nuclear-Electric Propulsion (NEP). 

The propulsion model can be modified to account for “limited” power sources 
such as Solar-Electric Propulsion (SEP) by constraining the power available to the engine 
using 


T) PA 

1 \ solar 1 0 1 array 


1 avail 


(7) 


where 

Pavail = electrical power available to the engines 

Po = reference solar incident power (power incident at 1 AU) in W/nr 
Aan-ay = solar array area in m 2 

T| S oiar = power converting efficiency of the solar array 

r au = vehicle distance from the sun (in AU) 

The power available can then be used to constrain the maximum thrust available at a 
given distance from the sun. 


F. NON-DIMENSIONALIZATION 

Non-dimensionalization (or scaling) is the reformulation of the problem such that 
all the optimizable parameters as well as the cost function assume values close to unity 
over most of the domain of interest. When not scaled properly, the problem may become 
“ill-conditioned” and the effects of numerical artifacts as common as round-off error can 
begin to impact the solution. Very poorly scaled problems may even behave as though 
there are singularities when in fact there are none. Problems that have not been scaled 
are susceptible to various ailments ranging from slow convergence to complete failure to 
converge. Betts [Ref. 5] gives many tips for the scaling of problems and points out that 
even when a problem is well scaled at one point, it may be scaled poorly at another point. 
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A simple non-dimensionalization method is to simply divide a state variable x by 
a representative reference value for that variable [Ref. 6]. The low thrust problem may 
be normalized in this way using a modified version of the approach implemented by 
Bryson [Ref. 1], Let a distance unit be defined as 1 AU and let a velocity unit be defined 
as the circular velocity of the earth. Thus we have 



where 

Udist = distance unit (1 AU) 

U V ei = velocity unit (circular velocity at 1 AU) 

Using this convention, a given radius r x could be non-dimensionalized as follows 



( 10 ) 


where an over-bar has been used to denote a non-dimensional variable. Continuing, the 
natural unit of time is simply 

U (ii) 

^ time T j V A A / 

U vel 

where 


Utime = time unit 

Note that for the above choice of distance and velocity units there are 2c time 
units per orbit. Now choosing the initial vehicle mass as our mass unit we have 

U mass = m 0 (I 2 ) 

where 


U m = mass unit 
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which completes the set of fundamental units. These particular choices of units are 
known collectively as canonical units. [Ref. 7] One advantage of using canonical units 
for astrodynamic problems is that the gravitational parameter p becomes one, 

simplifying many equations. A second advantage is that if only the above units (or 
combinations of the above units) are used to dimensionalize all the states, then the 
equations of motion take on exactly the same form as their dimensional counterparts 
provided there are no constant factors in any term which must be corrected (as is the case 
with the gravitational parameter which fortuitously is “corrected” to one). While it is true 
that the combinations of the above relations can be used to non-dimensionalize any 
parameter, there is no guarantee that the result will also be well scaled. In fact, the low 
thrust problem illustrates this fact nicely for particularly low thrusts. 

Consider a 600 kg vehicle with one NSTAR engine producing 92 mN of thrust. 
Using the units above, it would seem reasonable to assume that 


U 


thrust 


_ ^ mas dist 


(13) 


Substituting the values of the canonical we get 


U 


thrust 


(600kg) (1.496 x 10 11 m ) 

---- 5 -- = 3.9 N 

(5.023x10 6 )~ 


(14) 


Thus for this to be a useful unit for non-dimensionalizing thrust, we require our vehicle’s 
thrust output to be on the order of 4 N. For realistic low thrust systems, this value of 
thrust unit may be too high. A better choice would be to define a thrust unit as the 
maximum thrust produced by the engine. This guarantees that the normalized thrust 
magnitude will fall between zero and one but requires the equations of motion to be 
modified slightly as follows: 



(15) 


_ v , 

dt 7 


(16) 
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( 17 ) 


where 


dv r _ v, 2 1 + zT sinri 

dt 7 7“ in 

dv, _ _ v r v, zT cosrf 

dt 7 in 

dm _ zT 
dt V 

e 


^ thrush time 

dist ^d tel 


(18) 

(19) 


( 20 ) 


Clearly multiplying all the non-dimensional thrust terms in the equations of 
motion by this factor effectively accomplishes the same objective as normalizing thrust 
with purely canonical units. However, the difference is that domain of the thrust control 
is now well-scaled, and the conversions are taking place internal to the equations where 
proper scaling is not as important. 
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HI. AEROCAPTURE MODELING 


A. COORDINATE SYSTEMS 

The coordinate systems chosen for aerocapture are the same as those used in 
existing developed tools such as AC APS [Refs. 8,9] such that data could more easily be 
shared between programs. 

1. Position Vector 

The vehicle’s position vector is defined in the planet centered fixed (PCF) 
reference frame. This coordinate system is similar to a planet centered inertial (PCI) 
coordinate system except that the primary axis is aligned and rotates with a meridian line 
on the planet’s surface [Ref 7]. The angular velocity of the PCF frame with respect to the 
PCI frame is denoted by £2. As shown in Figure 2, the position vector is defined as 

r = [r,0 ,<f>] T where 

r = distance between origin and vehicle center of mass 
0 = longitude as measured positive in an easterly direction 

(f> = latitude measured positive up from the equator 
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Figure 2: Aerocapture Position Coordinate System 


Because the targeting of orbits with specific right ascensions is beyond the scope 
of this thesis and the equations of motion are invariant in longitude, it is convenient to 
define the longitude at atmospheric entry to be zero. 


2. Velocity Vector 

The velocity vector is measured in polar coordinates relative to the REN (Radius- 
East-North) frame. As shown in Figure 3, the REN frame is centered on the vehicle with 
the primary axis aligned with the extended radius vector. The tertiary axis is orthogonal 
and aligned toward the northern pole of the spherical coordinate system and the 
intermediate axis is mutually orthogonal and oriented easterly, completing the right- 
handed triad. 
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Radial 



Figure 3: Aerocapture Velocity Coordinate System (REN Frame) 

Within this frame, the velocity vector is represented by a magnitude and two 
angles, \ PCF - [v,\|/,y] T where 

v = vehicle speed 

\|/ = heading angle as measured counter-clockwise from easterly 

y = flight path angle as measures positive up from E-N plane 

It is important to note that the velocity vector represented in this frame is non- 
inertial unless Q, = 0 in which case PCF = PCI. For reasons that will become apparent 
later, it is sometimes convenient to express the vehicle’s inertial velocity in the local 

REN frame. This will be denoted by V PCT = [V, T, r] ' where capital Greek letters have 
replaced their corresponding lower-cased non-inertial symbols. 
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3. 


Frenet Frame 


The Frenet frame (sometimes referred to as the NTW frame [Ref 7]) is useful for 
resolving vectors into tangent, normal and bi-normal components with respect to the 
vehicle’s path. At a given time, the tangent component is tangent to the vehicle trajectory 
(alternately parallel to the velocity vector). The normal component lies in the plane of 
the orbit but perpendicular to the tangent component and the bi-normal component is 
normal to the orbit completing the right-hand triad. These components are identified by 
the subscripts s, n, and w respectively. The Frenet frame is a natural choice for resolving 
aerodynamic forces on a body as drag always acts in the negative direction tangent to the 
flight path and the lift vector will always lie in the normal/bi-normal plane. [Ref 7] That 
is to say that for a rotating atmosphere (fixed to the PCF frame) the Frenet axes are 
equivalent to wind axes. 


B. EQUATIONS OF MOTION 

For the generation of optimal trajectories, a three degree of freedom (DOF) model 
is sufficient to describe the spacecraft dynamics. The state of the vehicle is represented 

by the 7-dimensional vector x = [r,0 ,<[> ,v t|/ ,y]' consisting of radial position, longitude, 

latitude, speed, heading angle, and flight path angle. The equations of motion for a non¬ 
rotating atmosphere with a generic gravity model are given by [Ref. 10] as 


dr 

— = vsiny 
dt 

(21) 

<70 vcosy cost)/ 
dt rcoscf) 

(22) 

<7(|) vcosysint)/ 

dt r 

(23) 

dv 

~~ = a s +8 s 
dt 

(24) 

a w + 8 W -L cos y cost)/ tan(f) 
vcosy r 

(25) 
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( 26 ) 


dy _ a„+g„ [ vcosy 
dt v r 


where in the above equations a s , a n , and a w are the external accelerations resolved in 
the Frenet frame. 


1. Rotational Effects 

For the higher-fidelity case of a rotating atmosphere, Eqs.(9) through (11) must be 
modified to account for non-inertial centrifugal and Coriolis forces. 


dv 

dt 


= a s + g s +cf v 


ckf a + s v „ 

—- = — 1 --——cos y cost)/ tancf) + cf + co 

dt vcosy r ' 


dy a„+g„ vcosy . 

— = - — +- - + ct +CO 

dt v r 


(27) 

(28) 

(29) 


where the cf and co terms represent centrifugal and Corilolis contributions respectively as 
given by 


cf, =Q. r cos(|) (siny cos(|) -cosy sintf) sint)/ ) 


cL = - 


Qrr 

vcosy 


•sin (f) cos(|> cost)/ 


co =2f2(tany cos(|)sint)/ — sin cf)) 


Q?r 


cf =-cost]) (cosy cost]) + siny sintf) sint)/ ) 


co = 2f2cost|) cost)/ 


(30) 

(31) 

(32) 

(33) 

(34) 


and Q, is the planet’s angular rotation rate about its pole. 
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2. Aerodynamic Forces 

The aerodynamic forces of lift and drag are given [Ref. 11] as 


L=q d SC L {a) (35) 

D = q t SC B (a) (36) 

where Cl and Cd are the lift and drag coefficients and are in general a function of angle 
of attack (a ) and S is a reference area associated with the aerodynamic coefficients. 
The angle of attack is defined as the angle between vehicle’s body axis and the velocity 
vector as shown in Figure 4, 


Radial 



Figure 4: Definition of Angle of Attack and Flight Path Angle 


The dynamic pressure, q c \ is given by 


16 



( 37 ) 


V* =^p(f)v 2 

where p is the atmospheric density which is generally a function of altitude as given by 
the atmospheric model. 

3. External Accelerations 

The terms a s , a n , and a w represent the external accelerations resolved in the 
tangential, normal and bi- normal directions. These accelerations are the result of the 
aerodynamic forces and vary with bank angle, 8 which is defined as a body rotation about 
the vehicle’s velocity vector as show in Figure 5. 



Velocity 
(or tangent) 


Bi-normal 

Figure 5: Definition of Bank Angle 
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These accelerations can be expressed as function 

of the controls, mass and 

aerodynamic forces as 


D 

(38) 

a — - 

S 

m 


Lc os 8 

(39) 

a » = 

m 


L sinS 

(40) 

a =- 

W 

m 


Note that the total acceleration on the vehicle can 

be determined from these 

equations by using 


\\ a \\ = \l a s 2 + a n 2 + a w 2 

(41) 

and the total g-load becomes 


g - load = U 

(42) 

So 



C. GRAVITATIONAL MODEL 


For a simple inverse-square gravity model, the gravitational accelerations can be 

converted from their more familiar spherical form 

ft 

g=g r = - 2 

r 

(43) 

go =0 

(44) 

g* =° 

(45) 

to the Frenet frame using Eqn. (215) becoming 


g„ =-gcosy 

(46) 

g s = ~g siny 

(47) 

g w =o 
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(48) 



Increased fidelity gravitational models, such as those including zonal harmonics 
can be incorporated in a similar manner [Ref. 9] but is beyond the scope of this work. 

D. ATMOSPHERIC MODEL 

Anderson shows [Ref. 11] that if an atmosphere is assumed to be isothermal with 
constant temperature equal to the mean temperature, then atmospheric density varies 
exponentially with altitude. That is, 

(ft-jo) 

P( r ) = Po e H " (49) 

where 

p o = reference density 

H p = scale height 

h = altitude ( r-r p i anet ) 

ho = reference altitude 

The atmosphere is assumed to be fixed to the planet and rotates with the PCF 
coordinate system about the K PCF axis. The following table provides the above data for 
Mars and Neptune. 


Parameter 

Mars Value 

Neptune Value 

Po 

4.7x10 ^kg/m 3 

2.348x10 r 3 kg/m 3 

H P 

l.OxlO 4 m 

5.331xl0 4 m 

h 0 

4.9xl0 4 m 

0 m 

Atmosphere Limit 

12 5 km 

8 00 km 


Table 1: Atmospheric Parameters for Mars and Neptune 
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E. 


HEATING MODEL 


The thermal protection system (TPS) for a re-entry vehicle is typically sized to 
sustain both the heating rate and heat load (integral of heating rate) encountered by the 
spacecraft [Refs. 12, 13] These quantities take a maximum value at the stagnation point, 
the point on the surface of the vehicle where the freestream velocity is brought locally to 
zero. The heating rate at this point can be approximated by the Chapman equation 


<7 = 



(50) 


where 


r„ = vehicle nose radius 


C = stagnation point heating coefficient 

and N and M are constants for a given atmosphere associated with density and speed 
respectively. [Ref. 14] 


Parameter 

Mars [Ref. 15 ] 

Neptune [Ref. 17] 

C 

3.55x10 5 

7.9x10 5 

N 

0.5 

0.5 

M 

3.15 

3.0 


Table 2: Heating Rate Parameters for Mars and Neptune 


F. VEHICLE MODEL 

The vehicle data used was obtained from [Ref. 9] for consistencies with that 

author’s original work in aerocapture and is based loosely on the Mars Pathfinder shape. 

Due to the high cost to flight qualify new space flight components, most re-entry heat- 

shields have deviated very little since the early Viking designs. In light of this fact, the 

aerodynamic data below will likely be representative of possible future aerocapture 

vehicles. The aerodynamic coefficients listed below apply to Mars trajectories; lift and 
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drag coefficients at other planets should vary primarily due to variations in each planets’ 
atmospheric composition. Unfortunately aerodynamic data for a Neptune aerocapture 
mission is not available; therefore the values for Mars will be used with the 
understanding that these coefficients do not agree with the vehicle design. 

Following the atmospheric pass, the vehicle will require a small impulsive AV at 
apoapsis to place the vehicle in its final circular orbit. To accomplish this, a small 
conventional chemical engine using a bi-propellant comprised of Hydrazine (N 2 H 4 ) fuel 
and Nitrogen Tetroxide (N 2 O 4 ) as oxidizer is employed. This combination yields a 
specific impulse of approximately 330 seconds [Ref. 13] and provides sufficient thrust 
such that the AV may be considered impulsive. 


Parameter 

Symbol 

Value 

Vehicle Mass 

m 

568.5 kg 

Nose Radius* 

r n 

1 m 

Coefficient of Lift 

C L 

0.3024 

Coefficient of Drag 

C D 

1.6800 

Reference Area 

S 

5.52 nr 

Engine Specific Impulse 

Isp 

330 sec 

Engine Exhaust Velocity 

V e 

3234 m/s 


Table 3: Aerocapture Vehicle Parameters 

G. NON-DIMENSIONALIZATION 


The aerocapture equations of motion have been non-dimensionalized in the same 
manner used by [Ref. 9] in his original work in aerocapture and are repeated here for 
completeness. 


* Note: This value is not the actual measure of the nose radius; rather this value is chosen to be 
consistent with the stagnation point heating coefficient C in Table 3 
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_ r 

r 


U 


dist 


where U dist = r, net , the planet’s surface radius. 


U 


vel 


where U vel = __ ——, the circular velocity at the planet’s surface. 


planet 


( 51 ) 


(52) 


where U time = ^!£i and 


U, 


vel 


^time 


(53) 


_ m 
m — ■ 


U,„ 


(54) 


where U mKS =m 0 , the vehicle’s initial mass. Noting that the angles measured in radians 

are already conveniently non-dimensional and well scaled, Eqns. (21)-(23) may be 
presented in their non-dimensional form 


dr 

dt 


■ v siny 


dQ _v cosy cost)/ 
dt 7 c os(f> 

d§ _v cosy sint \f 
dt 7 


(55) 

(56) 

(57) 


The aerodynamic forces may be non-dimensionalize by first assuming that 

^ P 


U 


(58) 


density 
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where U density = p 0 , the scale density and using a non-dimensional reference area defined 
as 


5 = 


Noting that the dynamic pressure non-dimensionalizes as 


f U A 

mass 

. ^density^dist j 



the aerodynamic forces may be non-dimensionalized as 

L=q d SC L {a) 

D = q d SC D (oc) 

The external accelerations may now be written as 

D 

a s =- — 

m 

_ LcosS 

a n = -_ 

m 

_ LsinS 

a =- 

W — 

m 

The gravitational acceleration becomes 


where U gmv =— 

^planet 

Eqns. (27)-(29) become 


g = g, 



U 


grav 


dv _ _ 

~= = a s +g s +cf v 
dt 


(59) 

(60) 

(61) 

(62) 

(63) 

(64) 

(65) 

( 66 ) 


( 67 ) 
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( 68 ) 


chit a w + g w v , _ 

-—- - —cos y cos\)/ tan (j) + c/ + co 

dt vcosy r 


by 

dt 


a, + e„ vcosy _ 

+ —+ cf + CO 
v r 


(69) 


where the centrifugal and Coriolis forces have been non-dimensionalized using 

n=nu lime ( 70 ) 

such that 


£2 2 r cos(|) (siny cos(|) - cosy sintf) sint)/) 

(71) 

n 2 r 

c/ =-sin(f> costf) cost)/ 

(72) 

v cos y 


= 2£2(tanycos(|)sin\|/ — sin (f>) 

(73) 

Q. 2 T 

cos(f> (cosy cos([)+ siny sin([) sint)/) 

(74) 


v 
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IV. SOLVING OPTIMAL CONTROL PROBLEMS 


A. PRELIMINARIES 

An optimal control problem is the task of solving for the state and control 
histories of a system subject to constraints while minimizing (or maximizing) some 
performance index. Put mathematically (as presented by [Ref. 18] and repeated here for 
completeness), given a system with dynamic constraints such that 

x(x)=f(x,u,x;p) (75) 

where 

x = vector of states that completely describe the systems at any x 
u = vector of controls 

x = independent variable (usually but not necessarily time) 
p = vector of static parameters 
f = vector of dynamic equations 
subject to additional path constraints 

g / <g(x,u,x)<g„ (76) 

g = vector path constraint equations 

gi = vector of lower path bounds 
g u = vector of upper path bounds 
and boundary conditions (or point constraints) 

e , < e (x(x 0 ),x(x / ),x 0 ,x / )<e„ (77) 

e = vector path event condition equations 
ei = vector of lower event condition bounds 

e u = vector of upper event condition bounds 
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as well as bounds on states and controls 

x ; <x(x)<x u (78) 

u,<u(x)<u„ (79) 

in order to minimize a performance index of form 

x / 

J (x(»),u(»),x 0 ,x / ) = £'(x(x 0 ),x(x f ),x 0 ,x y ) + Jf(x(x),u(x),x )dz (80) 

to 

E = scalar cost function evaluated at the boundary times (event cost) 

F = scalar cost function evaluated over the entire time history (integral cost) 

When Fe0, the cost function is said to be in Mayer form. When Ee0, the cost 
function is said to be in Lagrange form. When E, F <t 0 the cost function is said to be in 

Bolza form. 

A mathematical construct known as the Hamiltonian may be created by adjoining 
a vector of Lagrange multipliers to the dynamic constraints and adding the integral cost 

H (x,u,A,,x;p) = F(x(x),u(x),x;p)+X T *f (x(x),u(x),x;p) (81) 

(Note: the Lagrange multipliers associated with the dynamics are also referred to as 
“costates”). The Pontryagin Minimum Principle (PMP) [Ref. 1] states that the optimal 
control history satisfies 

u = argmin H subject toucU (82) 

where u is the optimal control and U is the domain of u. 

This should be recognized as a static optimization problem at each instant of time 
for which the Hamiltonian itself is now the performance index minimized. We can now 
form an additional construct known as the augmented Hamiltonian (or Lagrangian of the 
Hamiltonian) defined as 

H 1 (p,A,,x,u,x;p) = H (x,u,?i,x;p)+p T -g (x(x),u(x),x;p) (83) 
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The Karush-Kuhn-Tucker conditions state that when evaluated at the extremal 


controls, the following are true 


du 

b t g = o 


<0 

& = 8 

II IV 
o o 

8 = g u 

&0 

VI 

so 

VI 

$0 

any 

3 

^0 

II 

^0 


(84) 

(85) 


( 86 ) 


where Eqns. (84) and (85) are referred to as the gradient normality and complementary 
conditions respectively. Eqn.(86) will be useful later for verifying the switching structure 
of the controls. 


B. SOLUTION METHODS 

Methods for solving optimal control problems can generally be separated into two 
groups, indirect and direct methods [Ref. 5 p. 85]. 

Indirect methods tend to produce greater accuracy and faster solution times. 
However, the problem formulation for an indirect method is considerably more complex 
as the user must provide additional information such as the equations for the costate 
dynamics as well as the gradient of the Hamiltonian with respect to the controls. 
Moreover, one must provide a reasonably accurate guess for the controls, states and 
costates. Unfortunately, the abstract nature of the costates makes this guess-work 
something of an art. 

By contrast, direct methods tend to be much more robust and can converge to an 
optimal solution even when seeded with a poor cr infeasible guess. Direct methods 
require less prior work on the part of the user as the user need only supply the dynamic, 
path and event constraint functions as well as the cost function. Direct methods work on 
the premise that a continuous problem may be approximated through careful 
discretization as a large number of point constraints, thus reducing the problem to a 
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single large Nonlinear Programming problem (NLP). Unfortunately, to increase the 
accuracy of the solution, greater resolution is required leading to larger problems that can 
be very slow to converge. 


C. DISCRETIZATION 

Approximation theory tells us that the worst (with regard to minimizing the L 
error norm) discretization scheme is to sample a function at equal intervals [Ref. 6]. 
Other methods of discretization include Hermite- Simpson and Sine methods [Ref. 20]. In 
fact, discretizing a function at the Legendre-Gauss-Labatto (LGL) points minimizes the 
Lo error norm for a given number of nodes. LGL points are characteristic in that they 
have a higher density distribution at the end points of a function and becoming sparser in 
the interior regions. By discretizing the problem in this manner, a direct solution may be 
obtained with either more accuracy or less computation time. 


D. DIDO 

DIDO is an application package [Ref. 18] for solving dynamic optimization 
problems in a friendly, easy-to-use MATLAB environment. DIDO employs a powerful 
direct Legendre pseudospectral method that exploits the sparsity pattern of the discrete 
Jacobian by way of the NLP solver SNOPT” [Ref. 21] and runs in both UNIX and PC 
environments. 

Problem formulation in DIDO is quite simple, with the user creating a set 
MATLAB functions to evaluate the problem dynamics, cost, path constraints and event 
conditions. These functions are tied together by an additional script which defines the 
upper and lower bounds for the states, controls, path constraints and event conditions as 
well as the initial guess. This guess need not be feasible; in most cases simply providing 
the estimated initial and final state value is sufficient. The simplicity of DIDO contributes 
to the rapid prototyping of solutions. An experienced DIDO user can generate optimal 
trajectories in days that previously would have required weeks or months. 
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The user defines the number of nodes (discretization points) for the problem and 
is afforded the option of “cold starting” or “bootstrapping” from a previous solution. 


E. COVECTOR MAPPING THEOREM 

Another advantage of the Legendre pseudospectral method is that is provides 
accurate costate and covector histories despite the fact that the adjoint equation are not 
supplied. [Refs. 22, 23] This power comes from the realization of the Covector Mapping 
Theorem (CMT). [Ref. 21] provides the following explanation of the CMT: 

The CMT may be articulated as follows: Given an optimal control 
problem, P , let P N denote its Legendre pseudospectral approximation 
where N is the order of the Legendre polynomial used in the 
approximation. Let P' denote the boundary value problem (BVP) arising 
from an application of the PMP to problem P , and P L denote the 
generalized root-finding problem obtained by applying the Karush-Kuhn- 
Tucker (KKT) conditions to problem P N . The CMT asserts that if P' is 
discretized by a Legendre pseudospectral method (i.e. an indirect method), 
then there exists an order-preserving map between these discretized 
covectors and the KKT multipliers associated with problem P NX . Hence, 
from the KKT multipliers, one can find covectors by the direct Legendre 
pseudospectral method as if one solved the problem by the indirect 
method. This is the CMT. 


The power of the CMT is that it provides the user with excellent tools for 
ascertaining the optimality of a given solution as we will see in the next section. 
Additionally, because the CMT provides accurate costate information, a DIDO solution 
can be uses as an extremely high quality guess for a more accurate indirect method. 


F. VERIFICATION OF OPTIMALITY 

When solving optimal control problems one is often challenged as to how one can 
prove optimality. [Ref. 21] discusses techniques for verification of optimal trajectories. 
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1. Feasibility 

Of primal concern is the feasibility of a solution. That is, does the solution satisfy 
the dynamic constrains of the problem? To ascertain this, the initial conditions are 
propagated using a Runge-Kutta algorithm (the MATLAB command ocle45 which 
implements the Dormand-Prince pair) using controls interpolated from the DIDO 
solution. If the DIDO state history and the propagated state history match (within some 
tolerance) feasibility is declared. 

2. Accuracy 

The next issue addresses the accuracy of the solution. Assuming the trajectory is 
feasible, the event conditions are evaluated using either the DIDO or propagated solution 
and examined to verify that all event constrains are satisfied (again, within a tolerance). 

3. Well-behaved 

Similarly, the solution is said to be well behaved if the path constraints (including 
state and control bounds) are obeyed over the entire trajectory. A solution that violates 
the path constraints may indeed be feasible; however it does not solve the problem at 
hand. 

4. Optimality 

Finally and most importantly, does the solution satisfy the necessary conditions 
for optimality? Recalling Eqn. (84), analytical expressions can be obtained by taking the 
partial of the augmented Hamiltonian with respect to each control. In some cases, these 
expressions can be rearranged into an explicit control law. Regardless, these expressions 
must hold true for there to be first-order optimality. These expressions can be verified 
through substitution of the states and CMT obtained costates and covectors. In the case 
where the necessary conditions can be re-arranged to solve for a control, these controls 
can be solved (referred to as CMT controls) and compared to DIDO solution controls. 
When the DIDO and CMT controls are in agreement, first-order optimality is declared. 

A second method for determining optimality is by examination of the 
Hamiltonian which is also constmcted by DIDO via the CMT. Recalling again Eqn.(84), 
the first integral of the necessary condition is required to be constant for Hamiltonians 
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that are not explicit functions of time. Thus the Hamiltonian can be plotted and its 
“flatness” used as a measure of optimality. 
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V. OPTIMAL LOW THRUST PROBLEM FORMULATION 


A. EVENT CONDITIONS 

For the orbit transfer problem, the vehicle must begin at Earth’s orbit distance 

r ( t o) = r Ear,h (87) 

and complete the trajectory at the target planet’s orbital distance 

r W = W (88) 

Assuming circular orbits with no ephemeris, there are an infinite number of 
equivalent solutions because the Hamiltonian is invariant with the vehicle’s angular 
displacement. Thus it is convenient to tie down the initial angular displacement to some 
value. 

e(*o)=0 (89) 


Assuming that the vehicle escapes Earth’s SOI with zero Q, the vehicle’s initial 
velocity must match that of Earth. Again, assuming circular planet orbits, the equation 
for a body’s radial and tangential velocities are 

v,(O = 0 (90) 



where p Sun = 1,327x 1O 20 is the Sun’s gravitational parameter. Moreover, initial 

vehicle mass must lie between the maximum capacity for the launch vehicle and some 
lower design bound. 

For the rendezvous case the vehicle’s final velocity vector must match that of the 
target planet thus 

E (t f ) = 0 (92) 
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(93) 






Sun 


' target planet 


For the case of a fly-by, the final velocity components would be left unconstrained. 

B. COST FUNCTIONS 
1. Minimum Time 

To minimize time, the performance index can be constructed in either Mayer or 
Lagrange form. In Mayer form, the cost function is simply 

J=t f (94) 

Formulated as a Lagrange cost, the performance index becomes 

J = { dr (95) 


2. Minimum Fuel 

For a single-staged vehicle, minimizing fuel is equivalent to maximizing the final 
vehicle mass. Maximizing a quantity is the same as minimizing the negative of that 
quantity, thus the Mayer performance index for minimum fuel is simply 

J = —m f (96) 


where m/ is the final vehicle mass. Alternately this performance index may be formulated 
as a Lagrange cost by recalling that the fuel flow rate for a constant specific impulse 
propulsion system is given by 


m 


fuel 



(97) 


so the performance index can also be expressed as 
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(98) 



dr 


C. CONTROLS 

The natural choices of controls for the low thrust problem are thrust magnitude 
and thrust angle. However, for some problems these choices lead to additional 
difficulties. Bounding the thrust angle as 0<r| < 2n can cause a jump discontinuity 
when the optimal control history passes through the boundary. Increasing the bounds to 
—2tc <r) < 2k yields multiple values of T| that are equivalent, again possibly contributing 
to numerical instability of the solver algorithm. 

To remedy this problem the equations of motion were altered such that the 
controls are the radial and transverse thrust components as well as thrust magnitude. 
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Thrust Mag 



Figure 6: Thrust Cone 


From Figure 6 clearly the thrust vector can be described by its radial and 
transverse components 

T r =T mag s inti (99) 

T , ~ T mag COST) 

and Eqns. (3),(4) and (5) can be rewritten with our new controls as 

fo, _ V M- , T r 

dt r r m 

fo, _ v , v , ! T , 

dt r m 


( 100 ) 

( 101 ) 

(102) 
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(103) 


Our control vector u is now 


dm 

dt 




where our controls are constrained within the following bounds 


-T <T <T 

max r n 

_T < j < j 

max t n 

o<r <t 

mag ma 


(104) 


(105) 

(106) 
(107) 


D. STATE BOUNDS 

Each state may be limited to a certain range of values over the trajectory, while 
others appear unbounded. As an example, the vehicle mass over the duration of the 
trajectory may not exceed its initial mass. Assuming that the vehicle’s only method of 
shedding mass is via thrust, then the dry mass of the spacecraft serves as a lower bound. 
Thus over the entire problem, mass is restricted to 

m dry ~ m ~ m initial (108) 

which might take a non-dimensional for of something like 

.1 <m<l (109) 

where the non-dimensional value for the lower bound depends on the vehicle’s dry and 
initial masses. 

Other states should be bound for more practical reasons. For example, a logical 
choice of bounds for heliocentric radius might be 

0<r<°° (110) 
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However, the dynamics [Eqns. (2)-(4)] become singular at r=0, thus this is a 
bad choice for lower bound. Similarly, although there is no physical reason to have an 
upper bound, placing a realistic upper bound of, say, 100 reduces the search space for the 
NLP solver. Thus better state bounds might be 

.1 < r < 100 (111) 

where the bounds have been expressed in canonical units. For an orbit transfer to Mars 
(at only 1.52 AU) this upper bound could be reduced reasonably further. 

Following similar reasoning, the remaining states take non-dimensional bounds of 


-10 <v,, <10 
-10 < v, < 10 
0<6 < 107t 

Thus the state bounds can be expressed in vector format as 


" .1 ~ 


r 


'100" 

0 
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10k 

-10 
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K 

< 

10 

-10 


K 


10 

.1 


m 


1 


( 112 ) 

(113) 

(114) 


(115) 


E. PATH CONSTRAINTS 

As currently bound, the control space incorrectly takes the shape of a rectangular 
solid. An additional constraint is needed to relate the three controls and correct the 
control domain. 

T, 2 +T, 2 -T mag 2 = 0 (116) 

In the control space u cz 9t 3 this path constraint represents the surface of a cone. 
Because the interior of the cone is not in the domain of u, the problem is not convex, thus 
leading to additional difficulties with attaining a rapidly converging solution. 
Interestingly, this problem can be convexified by restating the path constraint as 
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~T 2 <t 2 +T 1 _ T 2 <Q (H7) 

x mag — ± r 1 ± t ± mag — w V x x 7 / 

which represents a solid cone. Thus physics must be violated to convexify the problem. 
Fortunately, Pontryagin’s Minimum principle [Ref. 19] ensures that the optimal control 
lies on the boundary of the constraint surface. In essence, physics is intentionally 
violated to improve the convergence of the solution with the guarantee that the Minimum 
Principle will restore continuity of physical reality. 


F. NECESSARY CONDITIONS 

Recall from chapter IV that the Hamiltonian is defined 

H (x,u,A,,x;p) = F(x(x),u(x),x;p)+>i T »f (x(x),u(x),x;p) (118) 

Assuming any Mayer cost (that is F e 0), the Hamiltonian can be constructed as 


H=\(v r )+\ 


9 a,, 

V r ) 


v, [I T 

— -4+— 

r r m 


A ( 

+X 


_2Zi + ZlV?i ( T,m 


v r m J 


\ Ve J 


(119) 


The augmented Hamiltonian H ' is found by adjoining the path constraints to the 
above Hamiltonian as follows 


H f =K{v r )+\ 


fy, ^ 


\ r J 




V_iL + Zk 

r r 2 nr 


+ X„ 


V V 


T \ 


l + A,.. 


v r m J 


( t \ 
± mag 

\ Ve J 


+ ... 


( 120 ) 


[l 7 ,T +|i f + [i T + (J. ) 

r l r r r If t r l mag mag r rel \ r t mag / 


where , m., and u 7 are the Lagrange multiplies associated with the control bounds 

1 r 1 mag 

and p re/ is the Lagrange multiplier associated with path constraint relating the three 
controls [Eqn. (117)]. 

Applying the necessary condition for optimality by taking the partial derivative of 
the augmented Hamiltonian with respect to each control we obtain the following three 
necessary conditions: 

})H X 

^— = —+^+ 2^=0 ( 121 ) 
61 777 
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dH' ^ 

-^T = —+Mt +^re,T t =0 


(122) 


P)TT t ) 

—-= —^ + |i T -2u ,T =0 

-^y-. *T mag r^rel mag 

1 mag v e 


(123) 


These can be re-arranged to explicitly solve for the controls 


i (K 

T r =-^~ 

2 T,w m 


(124) 


1 

T < __ o71- 

2(1,, d m 


(125) 


T = 1 (-^ 2 l + m 

2,J v,^ 


(126) 


The normalized equations of motion had an additional factor relating the thrust 
units to the other canonical units. Thus in non-dimensional form Eqns. (124)-(126) 
become 


1 f z\ 

2 m 


(127) 


i r a,. 


2(M m 


r~+Mr f 


(128) 


T = 1 ( -+ M 

2u, v, + Mt - 


(129) 


where 


^ thrush time 
U dist U vel 


(130) 


These expressions will prove useful for verification of the DIDO solutions in the 
next chapter. 
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VL OPTIMAL LOW THRUST RESULTS 


A. PROCESS 

The code was first validated by solving the well-known minimum time low thrust 
transfer problem as presented in [Ref. 1], With this success in hand, the parameters were 
modified to better represent actual technology. It was at this time that the problems with 
the original non-dimensionaliztion scheme for thrust became apparent. Although 
minimum time trajectories were acceptable, feasible minimum fuel trajectories were 
difficult to obtain. The problem was reformulated using the improved non- 
dimensionalization scheme with better results. 

In general, low thrust problems took considerably less computation time than their 
aerocapture counterparts due to the fewer number of state variables and “slower” 
dynamics. For this reason, solutions were generally not bootstrapped, but instead solved 
completely each time from an initial guess. 


B. MINIMUM TIME RENDEZVOUS 

The minimum time, Mars rendezvous problem was solved for a vehicle powered 

by six NSTAR ion engines providing a total thrust capacity of 0.55 N. The initial vehicle 

mass was fixed at the maximum lift mass (C 3 of zero) for the Delta II 7326-9.5 which 

corresponds to 659.3 kg. If the initial mass is not fixsd at some value, initial runs showed 

that the solution would always choose the smallest initial mass available with the 

minimum propellant mass necessary to complete the trajectory, essentially minimizing 

the inertia the thruster must work against. As discussed in Chapter V, a zero C 3 departure 

trajectory can be modeled as beginning at the origin planet’s location, matching the 

planet’s tangential velocity and with zero radial velocity. Sixty nodes were used to solve 

the problem, requiring only 4.91 minutes on a Pentium IV PC (1.8 Ghz, 512 Mb RAM). 

This large number of nodes was more than sufficient to capture the details of the 

trajectory but the increased number of nodes was used to generate more aesthetically 

pleasing plots. The minimum time for transfer is 182.57 days during which time the 

vehicle consumes 270.2 kg of propellant for an arrival mass of 389.2 kg. Figure 7 shows 
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the non-dimensional state histories versus time. The symbols represent the DIDO 
solutions while the line represents the propagated solution (using the DIDO controls). 


Low Thrust States: Earth to Mars 



Figure 7: State History (Min Time to Mars Rendezvous) 


The radius vector smoothly and asymptotically transitions from 1 AU to 1.52 AU 
while the transfer angle increases almost linearly to 136.6 deg. Radial velocity increases 
to a maximum of 0.35 canonical units before symmetrically returning to zero to effect the 
rendezvous. Transverse velocity first slightly increases before beginning to decay 
matching Mars’ circular velocity of 0.81 canonical units. Mass of course decreases 
linearly due to the constant specific impulse and constant thrust profile as we shall see 
shortly. 

The heliocentric transfer orbit is shown in Figure 8 with the viewpoint from the 
solar ecliptic plane north pole. Again, the circles are the DIDO trajectory and the solid 
line is the propagated trajectory. The arrows are oriented with the thrust direction and 

scaled to the magnitude of the thrust vector. There is a noticeable switch approximately 
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half the distance between Earth and Mars where the thrust vector swings through from 
being predominantly outward to being predominantly inward. It is this change in the 
radial component of the velocity vector that begins to arrest the vehicles radial velocity to 
prepare it for a zero relative velocity arrival. 


Low Thrust: Earth to Mars 



Figure 8: Heliocentric Trajectory (Min Time to Mars Rendezvous) 

Figure 9 shows the control history during the trajectory. The DIDO controls are 
plotted as circles and the CMT derived controls are plotted as dots. The change in radial 
component of the velocity vector is readily apparent in the thrust angle control history 
just prior to the 100 day mark. As expected, the DIDO thrust magnitude for the 
minimum time trajectory employs maximum thrust for the entire trajectory. The DIDO 
controls and the CMT controls are in excellent agreement for the thrust anglealthough 
they do not coincide for the thrust magnitude. Despite the CMT controls telling us 
otherwise, the constant thrust solution is well known for the minimum time problem; thus 
it seems there may be an error in the way DIDO calculates costates and covectors or an 
artifact specific to this problem formulation. Figure 10 gives the Hamiltonian for the 
problem. Note that the variations in the Hamiltonian are on the order of 10'“, that is, the 
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Hamiltonian is very nearly constant. This fact, combined with the excellent agreement 
between the thrust angle DIDO and CMT controls suggest that the first order necessary 
conditions for optimality have been met. 


Control Histories 
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Figure 9: Control History (Min Time to Mars Rendezvous) 
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Figure 10: Hamiltonian (Min Time to Mars Rendezvous) 


C. MINIMUM FUEL RENDEZVOUS 

Next, the minimum fuel rendezvous was solved. The problem was set up 
identically to the minimum time problem with one notable exception, the initial mass 
event condition was free. As in the minimum time problem with initial mass free, the 
solution chose to depart Earth with the minimum amount of mass required to complete 
the journey. This is hardly in the spirit of planetary exploration, therefore the cost 
function was modified to instead maximize the final mass as in Eqn. (96). Thus in this 
case “minimum fuel” is a bit of a misnomer; instead, minimum fuel is implied by the true 
cost, maximized final mass. 

A second difficulty with minimum fuel trajectories is that they tend to gravitate 
toward extremely long trajectories in both time and path length. If both time and transfer 
angle are unconstrained, the vehicle will begin the trajectory with a very short duration 
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burn which slightly raises the aphelion of the trajectory. The vehicle then coasts for 
slightly more than a year before returning to perihelion where it again conducts a short 
bum, again slightly raising aphelion. This pattern is repeated over many years until the 
aphelion has finally reached the target planet’s orbital radius. These extremely long 
trajectories are difficult to capture using a direct method such as DIDO due to the 
comparatively few number of nodes typically used. To prevent this from occurring, the 
transfer angle was bound such that only type I orbits would be permissible (a type I orbit 
is an orbit with a total transfer angle between zero and 180 degrees, a type II orbit 
between 180 and 360 degrees, etc.) That is to say 

0<d(t f )<K (131) 

The optimal solution presented below is comprised of 60 nodes and 
required 7.85 minutes to converge. The total flight time was 0.695694 years (253.4 days) 
or 70 days longer than the minimum time trajectory. However only 119.1 kg of 
propellant are used giving an arrival mass of 540.3 kg. Note also that the final transfer 
angle is exactly 71 , the maximum allowable. Once again the propagated states (plotted as 
lines) are in excellent agreement with the DIDO states (plotted as various symbols). 
Thus the trajectory is at a minimum feasible. 
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Low Thrust States: Earth to Mars 



Figure 11: State History (Min Fuel to Mars Rendezvous) 


The heliocentric view of the trajectory [Figure 12] shows how this mass savings is 
attained. The vehicle begins by applying a near-tangential burn for approximately 40 
days before shutting off. This long bum places the vehicle in a transfer orbit whose 
aphelion exactly corresponds with Mars’ orbit radius. The vehicle coasts along this 
trajectory until around the 220 th day where it begins a shorter, sustained burn to increase 
the vehicle’s speed to match that of Mars. 
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Low Thrust: Earth to Mars 



Figure 12: Heliocentric Trajectory (Min Fuel to Mars Rendezvous) 

The next plot shows the control history throughout the trajectory. Again the CMT 
and DIDO controls match extremely well for the thrust angle but diverge for the thrust 
magnitude control. The DIDO thrust angle control (circles) can be seen to become more 
irregular in the middle portion of the trajectory when vehicle is not thrusting. This is 
expected as the thrust angle is not well defined without a thrust magnitude. Much more 
interestingly is the fact that the CMT thrust angle control history (dots) does appear to 
follow a well-defined, smooth curve even when the thrust magnitude is zero. Also note 
that although the thrust magnitude CMT and DIDO controls do not agree during the 
thrusting portion of the trajectory, they do agree during the non-thrusting portion. Thus it 
seems that the “switch” is only working in one direction. 
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Let us investigate this discrepancy further. Recall that the necessary condition for 
the thrust magnitude is 


mag 


2(i 


rel 


—— + (i r 

y mag 

V e J 


which can be rearranged as 


^-^mag^rel 


X. 


(132) 


(133) 


Because DIDO returns both the covectors with the solution (via the CMT) we can plug 
the DIDO solution states, controls and covectors into the right hand side of Eqn. (133) 
and check that the result matches the thrust magnitude covector as returned by DIDO. As 
we can see in the following plot, the two do not agree. 
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Figure 14: Comparison of DIDO and CMT Thrust Covector 
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We are left with two possibilities, either Eqn.(133) is incorrect, or the DIDO 
supplied thrust covector is incorrect. 

Despite the fact that we can not verify the DIDO controls by comparison 
to the CMT, we can still verify the thrust magnitude switches by plotting the switching 
function. Application of the KKT theorem tells us that 


\h. 


<0 
>0 if 
= 0 


- 0 


T =T 

mag max 

0 < T <T 

mag m 


(134) 


This implies that the controls switch between bounds at the zero crossings of the covector 
history. Figure 15 shows the thrust magnitude profile plotted above the switching 
function. The lower plot has been cropped to just either side of the x axis to better see the 
zero crossings. It is clear that the zero crossings of the lower plot coincide with the thrust 
switches as shown in the upper plot. 
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Figure 15: Switching Function (Min Fuel to Mars Rendezvous) 
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Finally, the Flamiltonian is plotted in Figure 16. Once again the Hamiltonian is 
approximately constant, verifying the necessary condition (via first integral). Note that in 
this case the value of the Hamiltonian is constant at zero. This is expected when time is 
not explicitly present in the cost function. 
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Figure 16: Hamiltonian (Min Fuel to Mars Rendezvous) 
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VII. OPTIMAL AEROCAPTURE PROPLEM FORMULATION 


A. EVENT CONDITIONS 

A spacecraft enters a planet’s SOI with a known velocity relative to the planet. 
This velocity, known as arrival V-infinity () determine the shape of the trajectory in 

the planet relative frame. When >0 , as in the case of aerocapture, this trajectory 

will be a hyperbola with the planet at the focus. For trajectories that do not intersect the 
atmosphere, the vehicle’s trajectory is Keplerian and easily solvable at any point. For an 
aerocapture trajectory, the vehicle will follow such a Keplerian trajectory from the arrival 
point to the upper limit of the atmosphere, sometimes referred to as the atmospheric 
interface. The atmospheric limit is defined as the altitude above the planet’s surface 
below which non-conservative forces such as atmospheric drag begin to perturb the 
trajectory. Thus only the atmospheric portion of the aerocapture trajectory need be 
optimized for two reasons: 

■ All portions of the trajectory outside the atmosphere have closed-form 
solutions that can be determined uniquely from the vehicle’s velocity 
vector evaluated at the atmospheric entry point (working backward to V 
infinity of arrival) and atmospheric exit point (working forward to the 
apoapsis). 

■ Above the atmospheric limit, the controls have no affect on the trajectory 
as there are no aerodynamic forces 
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Atmosphere Limit 
{Interface) 



Figure 17: Atmospheric Entry and Exit Points 


This being the case we can formulate our initial and final radius event conditions 
as 


( ^0 ) ^atm limit 

(135) 

| 

j 

ii 

(136) 


As was the case with initial angular displacement in the low thrust problem, the 
initial longitude is invariant. Thus it is convenient to define the initial longitude (as 
defined by atmospheric entry) to be zero. 

e(f„)=o (137) 
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The well know vis-viva equation [Ref. 2] defines the specific energy of a body at 
any point outside the atmosphere as a function of the vehicle’s inertial speed and radial 
distance from the central body. 


2 r 


However the arrival point is defined as r —> °o so 

vl 


’'arrival 


Due to conservation of energy, 


£ = £ 

arrival atm—in 


so substituting Eqn. (138) and (139) into Eqn. (140) we get 


v: v 2 . 

arrival atm—in 




r, 

atm 


where r atm is the atmospheric limit. This can be rearranged to 


V 2 -V 2 . I, 

~ atm— in , JA- 




- 0 


r „ 

atm 


(138) 


(139) 


(140) 


(141) 


(142) 


Eqn. (142) becomes the equality constraint relating the vehicle’s inertial velocity 
at atmospheric entry to V M . From the calculations in Appendix A, we can further 

relate the inertial velocity at atmospheric entry to the states as measured in the rotating 
frame using 

V 2 = v 0 2 + r 0 2 £2 2 cos 2 (f> 0 + 2r 0 v 0 i2cos(|) 0 cost)/ 0 covy 0 (143) 


where all state variables on the right-hand side are taken at atmospheric entry (initial 
time, also recall that lower-case velocity state symbols represent non -inertial 
components). 

Similarly, we require an event condition to ensure that the vehicle’s state at 
atmospheric exit is sufficient to target our desired apoapsis as show in Figure 18 below. 
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Figure 18: Targeting a Circular Orbit 


Again due to conservation of energy we have 


£ = £ 

atm—out apo 


(144) 


and using vis-viva we obtain 


y 2 y 2 

V atrn-o u t |X ' apo 

9 r 2 v 

^ 1 atm ^ 'apo 


(145) 


so V apo is needed as a function of the state vector at atmospheric exit. 


The angular momentum of a point mass is given by 

h =r xV = rFcos(r) (146) 

where T is the inertial flight path angle. Due to conservation of angular momentum 
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( 147 ) 



giving us an equality constraint relating the states at atmospheric exit to our desired 
apoapsis. Again, using the relations derived in Appendix A we can find V atm _ out and 

r atm _ ou , as functions of the state variables in the rotating frame using the following 
relations: 

V 2 = v f 2 +r 2 Ct cos 2 ^ + 2r / v / Qcos(|) / cost)/ / covy f (152) 

v / sin Y/ 

cos 2 y / (v^. 2 + f2 2 r / 2 ) + 2v / r / f2cosy / cost \r f cost^ 

where in this case all the state variables on the right-hand side are again taken at 
atmospheric exit. 

Note that these event conditions require that the apoapsis of the post-atmospheric 
trajectory exactly intersect the circular target orbit, also known as a Hohmann transfer. 
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Hohmann showed that a Hohmann transfer is the minimum-fuel transfer orbit between 
co-planar orbits whose ratio of radii are less than 11.94 [Ref. 7]. [Ref. 24] develops the 
event conditions necessary for a bi-elliptic transfer orbit where the apoapsis of the post- 
atmospheric trajectory has insufficient energy to reach the circular target orbit. Using the 
same approach, one could likewise develop the event conditions for the case where the 
post-atmospheric trajectory has excess energy and overshoots the target circular orbit. 
Unfortunately, these three cases (apoapsis undershoot, touching, and overshoot) can not 
easily be reconciled into one set of conditions because doing so would complicate the AV 
calculations in the next section, necessitating the introduction of an absolute value 
operator in the rocket equation Eqn. (156). Because the derivative of the absolute value 
operator is discontinuous when evaluated at zero, singularities are introduced in the 
Jacobian leading to serious numerical issues in solving the NLP. [Ref. 5]. For these 
reasons, only the “touching” case is considered in this work. 

The event conditions can be summarized in vector form as as 


0 

0 


< 




r (0 

e(0 


U; -U 

arrival 


atm—in ^ |*l 




r ( t f) 

cos(r ) 

V atm-out / 


-2p 


atm apo 


< 


0 

0 


(154) 


where V am-out and r a, m -„u, ^ given by Eqns. (152) and (153). 


B. COST FUNCTIONS 

1. Minimum Fuel to Circularize 

The simplest performance index for aerocapture is to minimize the fuel required to 
circularize the orbit. Since the event conditions require the post-atmospheric trajectory’s 
apoapsis exactly touch the target circular orbit, the delta-V can be obtained by subtracting 


the apoapsis velocity found using Eqn. (149) from the target circular orbit velocity 
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v„. 




( 155 ) 


such that 


AV=V. -V 

circ apo 


With the delta-V known, the rocket equation 


m pr OP =m 0 


f _A V } 

I, 


(156) 


(157) 


l — e 

v J 

can be used to find the propellant mass required to circularize the orbit. Thus our 
performance index can be expressed in Mayer form as 


J —m „ 


(158) 


2. Minimum Heat Load 

Knowing that there exists some relation between the total heat load and the 
required TPS mass, another reasonable performance index is to minimize the heat load. 
As one recalls, heat load is the integrated heating rate over the trajectory (Eqn. (50)) ; 
thus a natural choice is to use a Lagrange form cost function such as 

'/ 

J = Jg(r, v) dz (159) 

1 0 

However any Lagrange cost can be re-written as a Mayer cost through the 
addition of a state in the equations of motion. In this case, we could add heat load as a 
state, with dynamics given by Eqn. (50). 


3. Minimum Aerocapture Mass 

A more useful performance index would to minimize the total mass associated 
with aerocapture. The total aerocapture mass can be broken into three components, heat 
shield mass (also known as fore-shield mass), back-shell mass, and the propellant mass 
required to circularize the orbit. 
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prop 


( 160 ) 


The propellant mass to circularize can be found using Eqn. (157). The mass of 
the back-shell can be assumed to be relatively constant over a large range of heat. Thus a 
method is needed for determining the mass of the heat shield. 

A first-order approximation for heat shield mass can be made by linearly mapping 
heat load to heat shield mass. Thermal analysis of the heat loads expected to be 
encountered by Mars Smart Lander require an approximate heat shield mass of 40 kg for 
every 10,000 Joules of heat load [Ref. 25]. Using this single data point and assuming a 
20% margin, a crude mapping between heat load and heat shield mass can be expressed 
as 

tf 

m hea,shield = k \<1 ^ (161) 

1 o 

where k = 50k§ . 

10,000 J 

Since the back-shield mass is assumed to be constant, it can be neglected from the 
performance index. A minimum aerocapture mass performance index can now be written 
in Bolza form as 


‘f 

J = m prop + k\q dx (162) 

? 0 

It should be stressed that a much more accurate model could be obtained by 
mapping the heat shield mass to both heat load and peak heating rate. However, the 
addition of the peak heating rate term would require the cost function to be formulated as 
a Chebyshev minimax problem which is beyond the scope of this proof-of-concept. 
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c. 


CONTROLS 


For non-thrusting, fixed angle-of-attack aerocapture trajectories, the only control 
during the atmospheric pass is the bank angle of the spacecraft. As in the low thrust case, 
coding the control as an angle and bounding it as 0 < 8 < 271 leads to some difficulties 
due to the jump discontinuity between 0 and 2 k . As shown in later sections, the optimal 
bank angle history generally switches between 7t and 2rc (or zero). Opening the control 
bound to 0 <5 <371 allows for optimal bank angles that are entirely interior. However as 
in the low thrust case it was eventually realized that using sin6 and cosS as controls 
tended to yield faster and more consistently reliable solutions. Thus we have 



(163) 


where our controls are bound by 


-1 < sin 5 <1 

(164) 

-1 <cosS < 1 

(165) 


As in the low thrust case, this formulation of the controls requires the addition of 
a path constraint (given in the next section). 

An alternate control scheme incorporates the bank angle as a state variable of the 
vehicle and instead controls the bank angle. This allows for more realistic command 
response behavior. However, this scheme was not incorporated into this work. 


D. STATE BOUNDS 

Next suitable choices must be made to bound the states. Since only the 
atmospheric portion of the trajectory is of concern, there is no point in allowing the radius 
to be any larger than the atmospheric limit. Similarly, the vehicle has a hard lower 
bound at the planet’s surface. Considering this yields 


< r<r. 


(166) 


However, recalling Eqn. (49), the atmospheric model is only valid to the scale 
height below which the exponential atmosphere incorrectly decreases with decreasing 
altitude. Thus a better set of bounds is 
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( 167 ) 


planet 




atmos-limit 


Realizing that for a non-thrusting vehicle entering the atmosphere, the maximum 
speed will be near that at atmospheric interface. To account for the unlikely event of 
extremely steep entries the upper bound can be set to some arbitrary small multiple of the 
initial velocity. 


0< v<5v 


atmos -in 


(168) 


Note however that Eqns. (25), (26), (31), (33) have singularities at v = 0. This can be 
resoled by constraining the problem as 

8 <v<v„ (169) 

v — — atmos-in / 


where £ v is number slightly larger than zero. 

Note that equations Eqns. (22), (25), (31) also have singularities at <j) = ±K and 
y = ±Jt . We can prevent the singularities from being encountered by bounding these 
states as 


-K+% ^ ^ -£„, 

(170) 

—71 +£ y < y<7t -£ y 

(171) 


where again £ 0 and £ y are small numbers greater than zero. 


Finally consider longitude and heading angle. Since there are no singularities or 
physical bounds to contend with we are free to choose any convenient set of bounds 
sufficiently large to reduce the NLP search space. Note however that in both cases 


0 < 


0 

¥ 


< 271 


(172) 


is a poor choice. While we as humans intuitively know that these variables are 
continuous across these bounds, the dynamics are blind to this fact and instead see this as 
a discontinuity. This discontinuity can be resolved by simply opening the bounds 
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The bounds for the complete state vector can be summarized as 
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(174) 


E. PATH CONSTRAINTS 

The control space in U(z91 2 defined by Iqns. (164) and (165) is the area of a 
square. A path constraint is needed to correct this control space to that of a circle. This 
is accomplished using the trigonometric identity 

sin 2 8 +cos 2 S -1 = 0 (175) 

However this control surface leads to a convexity problem. Using the same 
technique employed for the low thrust control, the path constraint is modified to an 
inequality constraint 

-1 < sin 2 5+cos 2 5-1 < 0 (176) 

This transforms the control space to a circular disc of unit radius which is a 
convex surface. Again, this temporary violation of physics will be resolved by the 
Minimum Principle which forces the controls to the boundary of the control space, in this 
case, onto the unit circle. 

Other path constraints of potential importance to the mission designer are limiting 
maximum dynamic pressure, load-factor or heating-rate as given by Eqns. (37), (42) and 
(50) respectively. 
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F. 


NECESSARY CONDITIONS 

Assuming a Mayer cost (Fe 0), the Hamiltonian for the aerocapture problem is 


H = X r (vsiny) + L e 


f v cosy cost)/ ^ 


V 


+L 


+\ 


+K, 
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rcostf) 

A 




^vcosysint)/ ^ 


-gsiny + cf v 

m 


LsinS 


(177) 


v 


cosy cost)/ tan (|\> + cf +co 


mv cos y r 
LcosS g cosy vcosy 


-+ • 


mv 


+ cf y +co y 


Note that all of the aerocapture Lagrange costs are pure state costs meaning that they are 
not function of the controls. Since the partial of the augmented Hamiltonian will be taken 
in a moment, the assumption of a Mayer cost does not impact the formulation of the 
necessary conditions. 

The augmented Hamiltonian is 


H'=\ 


LsinS 
mv cosy 


+ \ 


LcosS 


v mv j 


+ |l sS sinS 


(178) 


+ 


(l c5 cosS +(i re/ (sin 2 S +cos 2 S -l) + pure state terms 


where all the pure-state terms have been grouped together and where (l sg , (l f5 and |l re/ 
are the Lagrange multipliers associated with the two controls and the path constraint 
relating them respectively [Eqn. (176)]. 

Applying the necessary condition for optimality by taking the partial derivative of 
the augmented Hamiltonian with respect to each control we obtain 


dH 1 
d (sin 8) 


mv cosy 


+ !L5 +2(i re /Sin5 =0 


(179) 


dH' _\L 
d (cos 8) mv 


+ (Ls +2(l re ,cos8 = 0 


(180) 


which can be re-arranged to solve for the controls 
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sinS =--(181) 

2mv\i rel cosy 2\i rel 

cos 8 =-—-(182) 

2mv[L re , 2ji reI 
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Yin. OPTIMAL AEROCAPTURE RESULTS 


A. PROCESS 

Problems were solved beginning with easier problem formulations before moving 
on to more challenging cost functions and higher fidelity. Solutions were obtained for 
non-rotating planets with fixed initial conditions and minimum fuel to circularize as the 
performance index. Next minimum heat-load trajectories were solved before attempting 
rotating atmospheres. With these successes in hand, the initial conditions were relaxed 
such that the velocity at atmospheric entry need only agree with the arrival V-infinity. 
This allowed for solutions with excess Vinfinity at arrival to be obtained. Finally, the 
solution for minimum total aerocapture mass was obtained. 


B. MINIMUM AEROCAPTURE MASS AT MARS WITH ZERO ARRIVAL 

V-INFINITY 

This case considered a vehicle arriving at a rotating Mars with zero excess 
velocity with minimum total aerocapture mass as the performance index. This initial 
arrival condition could occur in the case of a low thrust heliocentric trajectory whereby 
the interplanetary trajectory is a rendezvous. 

Although only 30 nodes were necessary to obtain a solution, using 50 nodes 
greatly increased the accuracy of the solution. The resulting performance index 
breakdown is given in Table 4. The total aerocapture mass (minus back-shell mass which 
is assumed constant) is only 19.14 kg. This low value can be attributed to poor modeling 
(poor choice for k ) in Eqn. (161) as well as the low heat loads experienced due to the 
zero V-infinity at arrival. 


Propellant mass 

10.63 kg 

Front-shield mass 

8.5 kg 

Total 

19.14 kg 


Table 4: Cost Function Breakdown (Zero Arrival V-inf) 
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Even taking these factors into account, the optimal aerocapture mass is of 
significant savings when compared to the required mass for a pure propulsive injection 
maneuver. Using the same main engine and bi-propellant modeled for the post- 
aerocapture circularization maneuver, the required propellant mass for the pure- 
propellant capture is 201.44 kg. Figure 19 shows the altitude, latitude and longitude 
histories during the atmospheric portion of the trajectory. The circles represent the DIDO 
solutions at the node points, whereas the line represents the Runga-Kutta propagated 
solution (note the strong agreement between the two indicating feasibility). As 
constrained, the trajectory begins and ends at the defined atmospheric interface of 125 
km. The minimum altitude of 70.14 km occurs at the 2.65 minute mark of the 12.24 
minute trajectory. Note that the initial latitude solution for the pass is very nearly zero 
(0.7 deg) and that the trajectory proceeds easterly. 
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Figure 19: Position History (Zero Arrival V-inf) 


Figure 20 shows the result of propagating the trajectory beyond the atmospheric 

interface. The spacecraft proceeds to an altitude of 295.4 km, only 4.7 km in error of the 
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targeted altitude. This error is attributed to the error in the final states at atmospheric and 
will be discussed more fully ahead. 


Propagated to Apoapsis 



Figure 20: Propagated to Apoapsis (Zero Arrival V-inf) 


The velocity state histories are given in Figure 21 and demonstrate the 
“backwards-S curve” in speed characteristic to aerocapture trajectories. The trajectory 
begins with an inertial speed of 4949.3 m/s and a relatively shallow flight path angle of 
-7.4 deg. The initial heading is only -0.09 deg, agreeing with the due easterly track in. 
Taken together, it is clear that the optimal solution is to fly in the direction of the rotating 
atmosphere at location of the atmosphere’s greatest velocity (the equator). This 
trajectory gives the minimum relative speed between the atmosphere and the vehicle, thus 
reducing the heating-rate which is proportional to V M . It is interesting to note that even 
when the problem is seeded with a westerly-tracking guess, DIDO still returns an easterly 
optimal solution. The inertial velocity at atmospheric exit is 3523.5 m/s for a total 
aerocapture delta-V of 1425.8 m/s. 
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Figure 21: Velocity History (Zero Arrival V-inf) 


The accuracy of the solution can be evaluated by comparing the DIDO terminal 
states with those from the propagator as shown in Table 5. The fact that the altitude 
indeed hits its target value (the only explicitly constrained final state) demonstrates the 
accuracy of the solution. The error that remains is likely due to these errors being less 
than the tolerances set in DIDO, thus causing the solver to exit. Furthermore, radius, not 
altitude is the actual state vector used in the solution. Redefining the equations of motion 
to use altitude as a state instead may allow for increased accuracy. 
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State 

DIDO value 

Propagator value 

Absolute Error 

Percent Error 

Hf 

125 km 

124.52 km 

-0.4762 km 

0.3824 % 

e / 

44.68 deg 

44.67 deg 

-0.0111 deg 

0.0248 % 

4 >/ 

0.11 deg 

0.11 deg 

0.0005 deg 

0.4335 % 

v f 

3274.5 m/s 

3272.7 m/s 

-1.8559 m/s 

0.0567 % 

¥/ 

0.09 deg 

0.09 deg 

0.0011 deg 

1.2560 % 

Y/ 

2.02 deg 

1.98 deg 

-0.0328 deg 

1.6523 % 


Table 5: Propagated Accuracy (Zero Arrival V-inf) 


Figure 22 shows the optimal bank angle history during the aerocapture pass. The 
blue circles represent the DIDO solution controls, whereas the red dots represent the 
CMT derived controls evaluated using Eqns. (181) and (182). For the bulk of the 
trajectory, the bank angle takes one of two approximate values, 0 deg (lift-up) or 180 deg 
(lift-down). This choice of extreme controls often occurs in optimal control solutions. At 
some points, the DIDO control solution actually oscillates between -180 deg and 180 deg 
which of course are equivalent. The same controls have been shifted to lie between 
0<5 <271 in Figure 23 to reduce this numerical irregularity. Again, note the excellent 
agreement between the DIDO solution controls and the CMT controls. The only 
significant divergence between the two occurs late in the problem when the vehicle is 
near its slowest and highest portion of the trajectory. In this energy state, the 
performance index is of greatly reduced sensitivity to the control (in other words, these 
errors are unimportant to the solution.) This agreement verifies that the first-order 
necessary conditions have been met strengthening the argument that the solution is at 
least a local minimum. Furthermore, the flatness of the Hamiltonian in Figure 24 proves 
verification of the first integral further contributing to the claim of optimality. 
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Figure 22: DIDO and CMT Controls (Zero Arrival V-inf) 
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Figure 23: DIDO and CMT Controls (0-2 %) (Zero Arrival V-inf) 
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The bank angle “switch” occurs at approximately 3.2 minutes, fully 30 seconds 
after minimum altitude passage. A possible explanation for this switching strategy is as 
follows: a steep, lift-up entry trajectory allows for deep atmospheric penetration with 
positive curvature. Once the “corner has been turned” at nadir and the velocity state is 
such that atmospheric exit is guaranteed, the vehicle flips lift-down essentially “clinging” 
to the atmosphere. This extends the time within the atmosphere but on the slower side of 
the minimum altitude where the heating rate penalty is less severe. 
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Hamiltonian vs. Time 
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Figure 24: Hamiltonian (Zero Arrival V-inf) 


The heating rate over the trajectory is given in Figure 25. The peak heating rate 
of 7.44 W/cm occurs at 140.18 sec (2.34 min). The total heat load was calculated as 

'y 

1700.38 J/cmr from the heating rate history using a trapezoidal integration scheme 
(MATLAB’s “ trapz” command). 
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Figure 25: Heating rate (Zero Arrival V-inf) 


The following figures show the time history of the body accelerations resolved in 
the Frenet frame and shows the total acceleration peaking at -0.85 g’s at time 158.94 sec 
(2.65 min), the vast majority of which is in the anti-tangential direction with small normal 
and bi-normal components. Not surprisingly, the peak dynamic pressure (shown in 
Figure 28) of 505.71 Pa occurs at the same time (recall from Eqn. (38) that the tangential 
acceleration is a function of drag and mass only). 
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Figure 26: Body Accelerations (Zero Arrival V-inf) 
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Figure 27: Total Acceleration (Zero Arrival V-inf) 
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Dynamic Pressure vs Time 



Figure 28: Dynamic Pressure (Zero Arrival V-inf) 


Finally, selected orbit elements are plotted versus time in Figure 29. One can see 
that the non-dimensional specific energy drops from zero energy (recall that the vehicle 
begins with zero arrival V-infinity) and that the eccentricity begins at a parabolic value of 
one before decaying to a near-circular value of 0.04 at atmospheric exit. The apoapsis 
begins at positive infinity (as it is singular for a parabola) and decays to a non- 
dimensional value corresponding to 300.0 km. 


76 




















Orbit Paramters (Non-Dimensional) 



Time (min) 

Figure 29: Selected Orbit Parameters (Zero Arrival V-inf) 


C. MINIMUM AEROCAPTURE MASS AT MARS WITH EXCESS 

ARRIVAL V-INFINITY 

This case considered a vehicle arriving at a rotating Mars with an excess velocity 
with minimum total aerocapture mass as the performance index. In this case, the initial 
arrival velocity is 5.2 km/s. This roughly equates to the arrival velocity expected from an 
impulsive Hohman transfer between Earth and Mars. 

The solution given below was obtained using the solution from Section B above 
and bootstrapping up by increments of 10 nodes until a satisfactory solution was obtained 
with 90 nodes. The cost breakdown for this solution is given in Table 6. Although this 
trajectory only requires one kilogram more of propellant, the required heat shield mass is 
slightly more than double the amount required for the zero arrival V-infinity solution. As 
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shall be seen, this is due to the significantly higher thermal loads placed on the vehicle 
during the atmospheric pass. 


Propellant mass 

11.81 kg 

Front-shield mass 

16.94 kg 

Total 

28.75 kg 


Table 6: Cost Function Breakdown (Excess Arrival V-inf) 


For this case, the cost savings when compared to a pure-propulsive capture are 
even more significant with the pure-propulsive capture requiring a whopping 386.4 kg. 
Figure 30 shows the altitude, latitude and longitude histories during the atmospheric 
portion of the trajectory. Note once again the strong agreement between the DIDO state 
history and the propagated solution. The minimum altitude of 58.5 km occurs 94.63 sec 
(1.58 min) into the 601.02 sec (10.02 min) trajectory. Again the optimal solution begins 
near zero latitude and tracks easterly with the rotating atmosphere. 
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Position History 



Figure 30: Position History (Excess Arrival V-inf) 


Figure 31 shows the continued propagation beyond the atmospheric exit point. In 
this case the propagated apoapsis reaches only 277.28 km, 22.72 km short of the targeted 
300 km altitude. This error is likely due altitude error at atmospheric interface (discussed 
shortly) magnified with time. 
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Propagated to Apoapsis 



Figure 31: Propagating to Apoapsis (Excess Arrival V-inf) 


Figure 32 gives the velocity state histories with an inertial entry velocity of 
7178.8 km/s and initial inertial flight path angle of -10.47 deg. The exit velocity is 
3517.1 km/s for a significant delta-V of 3661.8 m/s. 
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Figure 32: Velocity History (Excess Arrival V-inf) 


12 


As shown in Figure 33, this solution’s altitude trajectory begins with a steeper 
flight path angle and penetrates deeper into the atmosphere than the slower, zero arrival 
V-infinity solution. Note that although it penetrates deeper, the time of passage is 
significantly shorter. As we shall soon see, this contributes to a larger, but narrower 
heating-rate pulse which helps to reduce the total heat load. 
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Altitude vs. Time Comparison 



Figure 33: Altitude Profile Comparison (Excess Arrival V-inf) 


The accuracy and feasibility of the solution is shown in the following table. The 
error between the propagated solution and the DIDO solution is small enough to declare 
convergence. Note that although the final flight path angle has a large percentage error, 
the absolute error is less that 0.1 deg. Because the final value of the flight path is itself 
small, this makes the error appear large. 
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State 

DIDO value 

Propagator value 

Absolute Error 

Percent Error 

Hf 

125 km 

123.67 km 

-1.3323 km 

1.0773 % 

e / 

39.42 deg 

39.39 deg 

-0.0301 deg 

0.0763 % 

4 >/ 

0.20 deg 

0.20 deg 

0.0011 deg 

0.0763 % 

v f 

3268.11 m/s 

3262.32 m/s 

-5.7845 m/s 

0.1773 % 

¥/ 

-0.53 deg 

-0.53 deg 

0.0001 deg 

0.0270 % 

Y/ 

2.35 deg 

2.25 deg 

-0.0989 deg 

4.3985 % 


Table 7: Propagated Accuracy (Excess Arrival V-inf) 


The bank control history is given in Figure 34. The “bank-bang” nature of the 
control can be seen as the vehicle begins the trajectory lift-up before abruptly switching 
to lift-down at approximately 1.8 minutes. Again, the DIDO solution can be seen 
flipping back and forth between +180 and -180 which are of course equivalent. Both the 
switch as well as the agreement between the DIDO controls and the CMT controls is 
more prevalent in this case, likely do to the significantly larger number of nodes in the 
solution. Once again this satisfies the first order necessary conditions for optimality. 
Again the switch occurs after (but closer to) minimum altitude passage by approximately 
13 seconds. Hie optimality of the solution is further supported by the constancy of the 
first integral in Figure 35 (that is to say the flatness of the Hamiltonian to the 10" ). 
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Figure 34: DIDO and CMT Control Histoiy (Excess Arrival V-inf) 
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Figure 35: Hamiltonian (Excess Arrival V-inf) 


Hamiltonian vs. Time 
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The heating rate at the stagnation point is given in the following figure. The peak 

'y 

heating rate of 35.02 W/cnr occurs at time 72.76 sec (1.21 min) and the total heat load is 

'y 

3388.37 J/cm - almost double the value for the earlier non-rotating case. As mentioned 
previously, the duration of the heat pulse is quite short, less than a minute as measured at 
the half-maximum point. This contrasts with the non-rotating case where the half¬ 
maximum pulse width was approximately three minutes in duration. 
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Figure 36: Stagnation Point Heating Rate (Excess Arrival V-inf) 


The body accelerations are significantly stronger than those experienced in 
the previous case [Ref. Figure 37]. The tangential acceleration peaks at -4.5 g’s 
while the normal accelerations peak at 0.82 g’s. The total acceleration peaks at 
4.57 g’s at time 79.78 sec (1.33 min) [Figure 38]. Again, this corresponds exactly 
to the dynamic pressure peak of 2703.70 Pa [Figure 39]. 
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Body Accelerations 
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Figure 37: Body Accelerations (Excess Arrival V-inf) 
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Figure 38: Total Acceleration (Excess Arrival V-inf) 
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Dynamic Pressure vs Time 



Figure 39: Dynamic Pressure (Excess Arrival V-inf) 

The orbit parameters in the following figure are slightly more interesting 
than in the previous case. With this solution, the non-dimensional energy can be 
seen to begin well above zero, with capture occurring at 110.50 sec (1.84 min). 
The eccentricity decays from a hyperbolic 3.16 to a near-circular 0.04. Also the 
fact that apoapsis is undefined for parabolic trajectories is demonstrated as the 
apoapsis departs toward negative infinity to the left of the singularity before 
rapidly falling off from positive infinity immediately following the capture. 
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Orbit Paramters (Non-Dimensional) 



Figure 40: Selected Orbit Parameters (Excess Arrival V-inf) 


D. MINIMUM AEROCAPTURE MASS AT NEPTUNE WITH EXCESS 

ARRIVAU V-INFINITY 

The next case considered a minimum total aerocapture mass at Neptune with 
excess arrival V-infinity. The arrival V-infinity at Neptune was 9.42 km/s, significantly 
higher than that of Mars. A higher altitude of 1000 km was targeted for the final circular 
orbit and the atmospheric interface was defined as an altitude of 800 km. 

The problem was solved by first solving the zero excess arrival trajectory for a 
non-rotating atmosphere. This solution was used to bootstrap the excess arrival velocity 
which was in-turn used to bootstrap the case of a rotating atmosphere. 90 nodes were 
deemed sufficient to provide an accurate, optimal solution. Table 8 gives the cost 
function breakdown for this solution. The total required aerocapture mass is 328 kg of 
which 48.4 kg is propellant mass and 279.5 kg is heat shield mass. While significantly 
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higher than the aerocapture mass required at Mars, the total capture mass is 200 kg less 
than the 528.7 kg of propellant required for a purely propulsive capture maneuver. 
Again, the accuracy of the heat shield mass predictions is limited to the model used, 
which assumes heat shield mass scales linearly with heat load. 


Propellant mass 

48.4 kg 

Front-shield mass 

279.5 kg 

Total 

328.0 kg 


Table 8: Cost Function Breakdown (Neptune Excess V-inf) 


Figure 41 shows the position history of the spacecraft during the atmospheric 
portion of the trajectory. Note the excellent agreement between DIDO solution (circles) 
and the propagated solution (solid line). The trajectory begins and ends at 800 km of 
altitude (the defined atmospheric interface). The total pass requires 22.64 minutes with a 
minimum altitude of 256.5 km occurring 4.16 minutes into the trajectory. Again the 
solver chooses an initial latitude of nearly 0 degrees (equatorial) with an easterly track. 
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Figure 41: Position History (Neptune Excess V-inf) 

Figure 42 shows the result of propagating the DIDO solution beyond the 
atmospheric limit of the solution. The propagated apoapsis of 976.12 km is within 23.9 
km of the targeted 1000 km apoapsis altitude. 
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Propagated to Apoapsis 



Figure 42: Propagated to Apoapsis (Neptune Excess V-inf) 


Figure 43 gives the velocity state histories for the Neptune solution. The 
trajectory begins at a staggering inertial atmospheric entry speed of nearly 25000 m/s, a 
flight path angle of -9.35 degrees, and heading -0.1 degrees (due east). Inertial 
atmospheric exit occurs at 16144 m/s for an aerocapture delta-V of 8844 m/s. 
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Figure 43: Velocity History (Neptune Excess V-inf) 

The accuracy of the solution is presented in Table 9 which compares the DIDO 
and propagated terminal state values. Note that the “large” percentage errors in latitude 
and flight path angle are due to the small absolute value of the state. 


State 

DIDO value 

Propagator value 

Absolute Error 

Percent Error 

Hf 

800 km 

793.07 km 

-6.9257 km 

0.8733 % 

9/ 

48.32 deg 

48.31 deg 

-0.0116 deg 

0.0241 % 

<1>/ 

-0.00 deg 

-0.00 deg 

-0.0007 deg 

19.5726 % 

v f 

13736.80 m/s 

13732.23 m/s 

-4.5689 m/s 

0.0333 % 

V/ 

0.64 deg 

0.64 deg 

-0.0007 deg 

0.1067 % 

If 

1.51 deg 

1.46 deg 

-0.0495 deg 

3.3939 % 


Table 9: Propagated Accuracy (Neptune Excess V-inl) 
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The control history for the trajectory is presented in Figure 44. The DIDO 
solution for the first minute and last nine minutes appears somewhat erratic. 
Investigating this further, the covectors associated with the controls over these time 
intervals are of very small value, indicating a lack of sensitivity of the performance index 
to the bank angle in these regions. This corresponds with the physical explanation of 
reduced control effectiveness in the thinner upper limits of the atmosphere. In the 
thicker, lower atmosphere, the bank angle trajectory assumes the previously seen, lift- 
up/lift-down “bang-bang” type control with the switch occurring at approximately 5.5 
minutes, slightly after the minimum altitude point. The DIDO controls (circles) and the 
CMT derived controls (dots) are again in excellent agreement. 
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Figure 44: DIDO and CMT Control History (Neptune Excess V-inf) 


The thicker atmosphere and significantly atmospheric velocities contribute to 
extreme stagnation point heating for this trajectory. The maximum heating rate of 
264.59 W/cm occurs just prior to the minimum altitude point (Figure 45). The long time 
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duration of the pass results in a total heat load of 55,904 J/cm 2 driving up the mass 
requirements for the TPS. 
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Figure 45: Stagnation Point Heating Rate (Neptune Excess V-inf) 


Thermal protection is not the only engineering challenge presented by aerocapture 
at Neptune. The large aerodynamic forces lead to a peak total acceleration of 5.83 g 
(Figure 46). These high loads more akin to those experienced by fighter aircraft would 
require additional structural mass, further reducing available payload mass. This 
additional mass cost is not included in the modeling of this work. 
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Figure 46: Total Acceleration (Neptune Excess V-inf) 


E. MINIMUM AEROCAPTURE MASS AT MARS SUBJECT TO G-UIMITS 

The previous results at Neptune demonstrate the potential for high g-loads during 
aerocapture maneuvers. Rather than increasing the structural integrity of the spacecraft 
to survive the loads, another option is to simply constrain the g-limit to some smaller, 
more manageable value. The following solution shows partial results for the problem 
described by VIII.C (Minimum total aerocapture mass at Mars with excess arrival 
velocity) except that the tangential acceleration has been limited to 3 g. A 120 node 
solution was obtained by bootstrapping from the unconstrained solution results. The 
trajectories were similar with the constrained trajectory being slightly shallower (initial 
flight path angle of -9.73 degrees compared with -10.47) and longer in duration (total 
pass time of 11.4 minutes compared with 10.02 minutes). The position state histories are 
presented in Figure 47. The minimum altitude varies by only 3.5 km but occurs almost 
30 seconds later than the unconstrained case. 
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Figure 47: Position History (G-limited) 

The velocity histories are nearly identical, with a delta-V within 4 m/s of the 
unconstrained case. The only notable differences are some slight perturbations in the 
flight path angle and heading angle between times 1.5 minutes and 2.5 minutes during 
which time the spacecraft performs a bank maneuver to reduce the loads on the vehicle. 
A plot of bank angle versus time (Figure 49) shows this maneuver in greater detail. As 
before, the bank angle begins lift-up until approximately 1.4 minutes into the trajectory. 
At that time the vehicle banks instantaneously to lift-down for 30 seconds before 
switching back to lift-up for another 20 seconds. Finally, the bank angle switches lift- 
down for the remainder of the trajectory. 
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Figure 48: Velocity History (G-limited) 
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Figure 49: Control History (G-limited) 
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The first bank angle reversal is clearly reducing the tangential loads on the 
vehicle, as shown in the following plots of body accelerations. The first lift-down 
segment corresponds exactly to the time interval for which tangential acceleration is 
nearly constant at the constrained 3-g limit. This maneuver shows corresponding 
switches in the normal and bi-normal accelerations, coincident with the bank angle 
maneuvers. The second bank angle correction is more intriguing as it occurs after the 
tangential loads are decreasing in magnitude. This maneuver seems to be adding a slight 
increase in flight path angle to offset the period in which the flight path angle was 
relatively constant while lift-down (Figure 48). 
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Figure 50: Acceleration Histories (G-limited) 


The activation of the g-limit constraint can be further verified by examination of 
the covector associated with the constraint that is provided by DIDO. Figure 51 clearly 
shows that the constraint becomes active over the time period of the first bank reversal 
maneuver. 
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Figure 51: G-limit Constraint Covector History 


The heating rate is less than the unconstrained case, peaking at 30.39 J/cnr 
(unconstrained value of 35.02 J/cm 2 ) but the total heat load is higher at 3814.6 W/cm 2 
(unconstrained value of 3388.4 W/cm ). This is consistent with shallower trajectories 
which tend to have longer but smaller heat pulses. This higher heat load contributes to a 
slightly higher performance index as total heat load is mapped to heat shield mass in our 
model. The total aerocapture mass increases from 28.75 kg to 30.2 kg (Table 10). 


Propellant mass 

11.13 kg 

Front-shield mass 

19.07 kg 

Total 

30.20 kg 


Table 10: Cost Function Breakdown (G-limited) 


99 
























F. MAXIMUM AEROCAPTURE CORRIDOR SUBJECT TO HEATING 

RATE CONSTRAINT 

One of the difficulties with implementation of the aerocapture concept is the 
precision with which the spacecraft must guided to atmospheric interface. Aerocapture 
corridors are typically quite narrow, measuring only a couple degrees. The upper limit 
(that is to say shallowest angle) is normally defined as the shallowest initial flight path 
angle for which a lift-down bank angle profile will successfully meet the terminal 
conditions [Ref 9]. Similarly the lower limit is defined as the initial flight path angle for 
which a lift-up bank angle profile will meet the terminal conditions. Some preliminary 
work was done in investigating whether an optimal bank angle strategy can increase the 
aerocapture corridor width. 

The maximum aerocapture corridor problem was solved by separately solving two 
related problems: maximum initial flight path angle and minimum initial flight path 
angle; respectively: 


J=~ Yo 

(183) 

J= Yo 

(184) 


The difference between the minimum and maximum initial flight path angle is then the 
maximum corridor width. In addition, both solutions were subject to heating-rate 
constraints. Note that it is important that the initial conditions for each problem be the 
same so that the resultant trajectories can be fairly compared. In fact, if this is not done, 
the optimal solution for the two problems differ in initial heading by 180 degrees! 
Instead the minimum initial flight path angle solution was generated first, and its initial 
condition was used to constrain the initial maximum flight path angle case. Moreover, 
with no constraint on the final orbit inclination, the optimal solutions placed the vehicle 
in an equatorial orbit - not very desirable from a scientific point of view. However, for 
the equatorial case, the corridor-defining optimal bank profiles were constant lift-up or 
constant lift-down as per the assumption in the definition. This turned out to not be the 
case for trajectories with final inclinations other than zero. 
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Additional solutions were generated with a constraint on the final orbit inclination 
of the vehicle. Recalling that upper-case Greek letters represent the components of the 
velocity vector resolved in the inertial frame, the inclination of the vehicle is related to 
the aircraft states by: 

cos i = cos f cos 'R f (185) 

'y 

The maximum heat rate was set at 50 W/cmr and the targeted inclination was arbitrarily 
chosen to be 70 degrees such that 

cos(70deg) ^cosd^ cos'f / / <cos(70deg) (186) 

Relevant parameters for the two trajectories that bound the maximum corridor are 
presented in Table 11. 
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Minimum Initial FPA 

Maximum Initial FPA 

Initial Flight Path Angle 

-10.83 deg 

-8.66 deg 

Initial V-infinity 

5200 m/s 

5200 m/s 

Initial Latitude 

17.6 deg 

17.6 deg 

Initial Heading 

-73.4 deg 

-73.46 deg 

Delta-V 

3743.5 m/s 

3662.4 m/s 

Total Pass Time 

5.87 min 

11.82 min 

Minimum Altitude 

56.32 km 

66.52 km 

Max Dynamic Pressure 

3298.44 Pa 

1354.9 Pa 

Max Total Acceleration 

5.57 g 

2.29 g 

Max Heating Rate 

39.76 W/cm 2 

27.67 W/cm 2 

Heat Load 

3138.83 J/c 2 

4612.22 J/cnr 

Front-shield Mass 

15.69 kg 

23.06 kg 

Post-Aerocapture Propellant 

26.41 kg 

11.91 kg 

Total Aerocapture Mass 

42.1 kg 

34.97 kg 


Table 11: Max Corridor Boundary Comparison 


Figure 52 is a plot of the position state histories for the two trajectories (circles 
represent the maximum flight path angle entry while the plus symbols represent the 
minimum flight path angle entry). The trajectories begin at atmospheric interface at the 
optimal latitude of 17.6 degrees. The minimum flight path angle trajectory is steep 
(-10.83 deg), with only 5.87 minutes of total pass time compared with 11.82 minutes for 
the shallow maximum flight path angle case (-8.66 deg). The steeper trajectory naturally 
penetrates deeper into the atmosphere to an altitude of 56.32 km as compared with 66.52 
for the shallower trajectory. 
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Position History 



Figure 52: Max Corridor Position Histories 


The velocity histories are given in Figure 53. The total delta-V for the two 
trajectories are nearly equal at about 3703 m/s +/- 40 m/s. 


Velocity History 



Figure 53: Max Corridor Velocity Histories 
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The control histories for the trajectories are given in the following figure. CMT 
control histories were plotted in lieu of the DIDO solutions as they better demonstrated 
the arctangent characteristic of the bank angle schedule. Clearly the optimal bank angle 
profiles are not simply lift-up or lift-down as was the case for equatorial target orbits. 
The steep flight path begins lift-up and modulates to approximately 60 degrees during the 
ascent. In a like manner, the shallower entry begins lift-down and modulates to 
approximately 94 deg during the ascent. In neither case does the lift vector cross the 
horizontal plane; instead the lift vector seems to be used to control depth of penetration 
prior to the minimum altitude before switching to control heading (and hence target 
inclination) for the ascent. 
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Figure 54: Max Corridor Control Histories 


As expected, the peak heating rate for the steep entry is significantly higher than 
the shallow entry (39.76 W/cm" as compared to 27.67 W/cnr) although the shallower 
trajectory has the higher heat load (4,612 J/cm 2 compared with 3,138 J/cm 2 ) (Figure 55). 
Note that the peak heating rate was below the constraint value of 50 W/cm" 
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Figure 55: Max Corridor Heating Rates 


3500 


3000 


2500 


|2000 


= 1500 


Dynamic Pressure vs Time 


IV 

++ 

ax Dynamic Pr 

jssure = 3298.4 

1 Pa 


o Max y Q 
+ Min y 0 


+■ 4- 

4* 4 






4- 

4 






4 

4 

4 






4 ' 

o'* 

0 

4 o 

Max Dynam 

P°o 

4 - ° 

c Pressure = 1i 

54.90 Pa 



O 

* 0 

4 0 

4 , 0 

* 

+ °o„ 





4- u 
+? 

J 


f, 

°o Q 

O °°ooo f 

i0 °000000oo, 

ooooooaaaese 

'jxnxjxxsmsBm 


1000 


10 


12 


Time (min) 


Figure 56: Max Corridor Dynamic Pressures 
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The peak dynamic pressure (Figure 56) of the steep trajectory is more than double 
that of the shallow trajectory (3,298 Pa compared with 1,355 Pa) which leads to a similar 
disparity in peak total accelerations (Figure 57). The steep entry encounters a crushing 
5.57 g peak acceleration whereas the shallow entry peaks at a more manageable 2.29 g. 
Figure 58 resolves the aerodynamic accelerations into flight-path related components. 


Acceleration Mag vs Time 



Figure 57: Max Corridor Total Accelerations 
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Figure 58: Max Corridor Acceleration Components 


The significant differences in characteristics of the trajectories which define the 
boundaries of the maximum aerocapture corridor illustrate the difficulties imposed upon 
the design team. To utilize the entire available corridor, the vehicle’s structure must be 
sized to the more dynamic, steep boundary trajectory while the TPS mass must be sized 
for the larger heat loads of the shallow entry. Additionally, the TPS material selection 
will be dependent on the maximum heating rate of the steep trajectory. 
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IX. FORMULATING THE COMBINED LOW THRUST AND 
AEROCAPTURE PROBLEM 


Previous chapters detail the problem formulation for both low thrust and 
aerocapture trajectories taken separately. The attention is now turned to the problem of 
solving for the low thrust trajectory that terminates with an aerocapture maneuver. One 
approach for generating feasible integrated trajectories is to optimize the low thrust 
trajectory using its final state to derive the initial state for an optimal aerocapture pass. 
However, this formulation is not optimizing the problem from end to end. To find the 
true integrated optimal solution, we must simultaneously solve both trajectories. 


A. JUNCTION CONDITIONS 


Recall that a direct, discreet solution method reduces an optimal control problem 
to a series of dynamic constraints sampled at various nodes (times). Thus, the state 
history for a low thrust problem can be expressed as a matrix such as 


r (0 

r (0 

• r ( t j) 

e(f 0 ) 

0(0 

• 

V, (to ) 

V r(h) ■' 

•• v X t j) 

v ,('o) 

V r (0 •' 

■ v Xf,) 

m(t 0 ) 

m(/j) •• 

•• 


(187) 


Where the columns represent the values of the states at each node time and were j 
is the index corresponding to the final node of the low thrust solution. Similarly for 
aerocapture 



r (0 

••• r(t k ) 

e (0 

e (0 

- 0 ( 4 ) 

<K*o) 

00 

••• <K* t ) 


v (0 

••• v(t k ) 

¥ (0 

¥ 0 ) 

... \|/(r t ) 

y (0 

y (0 

Y ( 4 ). 
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(188) 






where k represents the final node index of the trajectory. 


If these two state history matrices were of the same dimension, they could be 
augmented forming a combined state history for the entire problem. However, note that 
the low thrust state history matrix has one less row (state) than the aerocapture matrix. 
This can be resolved by simply adding a “dummy variable” to the last row of the matrix. 


r ( f o) 

r (0 

• r ( f j) 

e(0 

e(0 

■ 

V r (*0 ) 

V r(h) •' 

•• v X t j ) 

v A) 

v Xh) 

■ v X r ,) 

m (t 0 ) 

m(/j) •• 

■ m ( 1 ,) 

X 

X 

• X 


(189) 


Our combined state history matrix now takes the desired form of a 6 by §+k ) 
matrix with the final column of the low thrust trajectory occurring at index j and the first 
column of the aerocapture trajectory beginning with column /+/. 


r ( r o) 

r (h) •• 

• Xt,) 

r (Vi) ' 

•• r (‘,J 

e (0 

e (0 •• 

• 0 h) 

9 (»/.,) ■ 

- »(>,.„) 

V r (t 0 ) 

v Xh) •• 

• V r(tj) 

♦(',. 1 ) ■ 

" 0 ) 

v,(t 0 ) 

v r(0 •• 

■ v X t j) 


" v {'m) 

m (t 0 ) 

m (h) ■■ 

• m (o) 

■ 

" v(v*) 

X 

X 

X 

v(v.) ■ 

•• tM 


However, we now must somehow match the physical meaning of the end of the 
low thrust trajectory with that of the beginning of the aerocapture trajectory. As 
formulated, the states take on very different meanings in the two portions of the problem. 
For example, the 5 h row of the low thrust portion represents the vehicle’s mass while in 
the aerocapture portion that same row represents heading angle. Thus we need some 
junction conditions to relate the variables on the two sides of the problem. 
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First, recall from previous chapters that the initial final radius of the low thrust 
portion and the initial radius of the aerocapture problem are fixed. That is 


r( tj ) = r 




planet 


T 

atmos limit 


(191) 

(192) 


Final velocities for the low thrust portion are left unconstrained to allow for non¬ 
rendezvous (excess arrival V-infinity) trajectories to be generated. The vehicle mass at 
the end of the low thrust portion remains subject to the physical limitation that 


m 


dry 


<m(t)<m(t 0 ) 


(193) 


Although mass is not a state for the aerocapture segment, its value is needed for 
the dynamics calculations. Assuming that there is no post low thrust staging this yields 


m , = m , 

aerocapture \ j 




(194) 


Since vehicle mass at arrival is on the same order as the initial vehicle mass at 
launch, there conveniently is no need rescale the problem as the normalizing units of 
mass may be chosen to be the same for both trajectory segments. 

Recall that for convenience the initial longitude of the aerocapture trajectory was 
set to zero so 


0 (%)=o 


(195) 


Similarly the initial velocity states of the aerocapture problem can be related in 
some manner to the terminal conditions of the low thrust problem. 


V 7+1 = f (' (tj ) > V, (tj ) , r (t j+1 ) ,<> (t j+1 ) ,y (t j+l ) ) 

(196) 

V /+ I =f(v r (tj ). V I (tj ), r (t j+ 1 ) ,<l> (tj + ,) ,y ( t j+l )) 

(197) 

=f( V r(tj)’ V '(tj)’ r (tj + 1 ) ^ ( V. ) • Y ( Vl ) ) 

(198) 
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To see exactly how these are related, the V-infinity of arrival for the aerocapture 
problem must be found. The heliocentric inertial velocity of the vehicle is give by the 
vector equation 


'vehicle planet f j ^vehicle f 

sun v J sun 1 J 


planet 


(199) 


which can be rearranged to find the inertial velocity of the vehicle relative to the planer 


{f 


1 =\V ) 

I vehicle 


vehicle ( . 

J planet 


| ^ planet } 


( 200 ) 


With respect to the variables used in the low thrust portion, the velocity vector of 
the vehicle with respect to the sun is 



and the circular velocity of the planet is given by 



0 


IT, 


planet 


( 201 ) 


( 202 ) 


Thus the velocity of the vehicle relative to the planet can be expressed as 


1 v vehicle J , 
v J planet 


V r 

v t - 

I 1 sun 



\ V 

y planet 


(203) 


The magnitude of this vector is 





f 

-3 

1 y vehicle I , 
v J planet 

= V„ = t 

arrival 4 

lv,. 2 + 

v t 2 - 

M' sun 

\ v 
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\ 

planet 


(204) 


The square root term leads to possibility of a singularity so it is more convenient to use 
the square of the V-infinity at arrival given as 




( 


+ 



V 



(205) 
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At this point, the inertial arrival velocity magnitude in the planet frame at the 
beginning of the aerocapture trajectory is known. Given the small size of the target 
planet relative to the scale of the heliocentric trajectory, small course corrections far from 
the target planet allow for any point on the B-plane to be targeted (the B-plane is a 
reference plane used for interplanetary targeting). This means that the initial latitude, 
longitude, velocity, heading and flight path angle sates for the aerocapture problem is are 
essentially free, provided they all take values that are consistent with the arrival V- 
infinity. 

Using vis-viva, the magnitude of the velocity vector at the atmospheric interface 
(where we begin the aerocapture trajectory) is related to the arrival V-infinity. 

-f^ + JL = o (206) 

2 r 


Again referring to Appendix A, the velocity components in a rotating REN frame 
may be related to the arrival V-infinity, thus we can now relate the V-infinity as 
calculated from the final states of the low thrust trajectory to the initial states of the 
aerocapture trajectory as 

fc)|f =v(t, + i) 2 +r(t, + i) 2 Q 2 C os 2 ^)(i, +1 )--- (2Q7) 

+ 2r (t j+l ) v (t j+1 ) f2cos(^ (t j+l ) cosy (t j+l ) cosy (t j+1 ) 


Substituting Eqn. (205) into the above equation yields the following important junction 
condition 


(U 


+ 




M 


A) 




(208) 


+ 


2r (t j+1 ) V (t j+l ) a costf) (t j+l ) cos y (t j+l ) cosy (t j+l ) 


However, recall that not only were the states different between the two 
formulations, but they were scaled markedly different as well. This can be resolved 
through the addition of a constant conversion factor 
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(‘,f 


+ 


(t ) 2 _ 

u ~H 


) 


^ dist _ lowthrust 
y ^dist _ aerocapture j 


v (Vi) >r (Vi)^ 2 cos >(Vi)- 


(209) 


+ 2 r (t j+1 ) v (t j+1 ) D. cos(^ (t j+l ) cosy (t J+ ,) cos y (t J+] ) 


Eqns. (191)-( 193), (195), and (209) complete the junction conditions for the 
combined problem. 

B. COST FUNCTIONS 

1. Minimum Total Propellant 

A relatively simple to implement performance index for the combined problem is 
to minimize the total propellant required for the mission. This can be accomplished 
despite the propellants being of different types with different I sp . This combined 
propellant cost can be expressed in Mayer form as 

J = — (m . , +m . ) (210) 

V proplowthrust prop circ / v 7 

where m proplowthrust is the propellant mass consumed during the low thrust portion of the 
trajectory as given by 

^proplowthrust "•(<„) (211) 

and m propcirc is the mass of propellant needed to circularize the post-aerocapture transfer 
orbit. 


2. Maximum Scientific Mass 

The total maximum scientific mass at arrival can found by taking the total 
propellant cost above and adding to it the estimated mass of the heat shield. Using the 
same mapping between heat load and front-shield mass used in chapters VII and VIII this 
can be expressed as a Bolza cost as follows 
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0 proplowthrust propcirc 


l k 

) + £ J g t/r 


( 212 ) 
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X. SOLUTIONS TO THE COM BINED LOW THRUST AND 
AEROCAPTURE PROBLEM 


A. MINIMUM TOTAL FUEL SOLUTION 

The minimum total fuel solution was obtained using 40 nodes of resolution for 
each segment of the trajectory. Obtaining feasible trajectories proved rather difficult. 
Increasing the thrust capacity of the low thrust vehicle helped to obtain solutions at the 
expense of realism. In addition, the rotation of the target planet was set to zero to simplify 
the solution and reduce computation time. This was accomplished by simply assuming 
that the vehicle employed 10 NSTAR engines. As show in Figure 59 (where again the 
circles are the DIDO solution and the solid line is the propagated trajectory) the trajectory 
begins with zero Q and a 700 kg vehicle and arrives at 1.52 AU 291 days after launch. 
In the process, 61.58 kg of propellant are consumed for an arrival mass at Mars of 
638.42 kg. The close agreement between the DIDO solution and the propagated solution 
verify the feasibility of the solution. 
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Low Thrust States: Earth to Mars 



Figure 59: Low Thrust State Histories (Min Total Fuel) 


The heliocentric trajectory is shown in Figure 60. An initial burn of 
approximately 48 days increases the semi-major axis of the transfer orbit until the 
aphelion intersects Mars’s orbit. The small deviations between the DIDO and propagated 
solutions are more visible in this presentation. 
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Low Thrust: Earth to Mars 



Figure 60: Heliocentric Low Thrust Trajectory (Min Total Fuel) 


The vehicle arrives at Mars with an inertial arrival V-infinity of 2626.2 m/s and a 
relatively shallow flight path angle of -7.62 degrees. Note that this is a relatively shallow 
flight path angle but the correspondingly small arrival V infinity ensures that the vehicle 
will not “skip” off the atmosphere. Figure 61 and Figure 62 depict the vehicles states for 
the aerocapture portion of the trajectory. A minimum altitude of 70 km is reached 3 
minutes into the 13.65 minute trajectory. The pass yields a total delta-V of 2078 m/s. 
Again the DIDO solution and the propagated solution are in excellent agreement. The 
transfer orbit is highly circular with an eccentricity of only 0.03. This contributes to the 
low propellant mass required to circularize of only 12.1 kg. 
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Position History 
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Figure 61: Aerocapture Position States (Min Total Fuel) 
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Figure 62: Aerocapture Velocity States (Min Total Fuel) 
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The bank angle history is given in Figure 63. Unfortunately, current versions of 
DIDO do not return the covectors for problems with internal knots such as this one 
making verification difficult. For this reason, the CMT controls can not be shown for 
comparison. Similarly the Hamiltonian is not available for verification of the first 
integral. 
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Figure 63: Aerocapture Control History (Min Total Fuel) 


Figure 64 gives the heating rate over the trajectory. The peak heating rate of 12.2 

9 9 

W/cnr occurs 2.14 minutes into the trajectory and the total heat load is 3054.3 J/cnr. 
Using 50 kg of TPS mass per 10,000 J/cm 2 of heat load, this corresponds to 
approximately 15.3 kg of front-shield mass. 
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Figure 64: Aerocapture Heating Rate (Min Total Fuel) 


Unfortunately, DIDO does not currently return covectors for problems with 
interior knots such as this mixed-dynamics problem. Without the co vectors or 
Hamiltonian, verification of the optimality of the solution is more difficult. In principle, 
the interior event conditions of the mixed-dynamic problem could be used to formulate 
two optimization problems, essentially breaking the problem back into its parts. Each 
optimal control problem could then be solved, and the covectors and Hamiltonians 
exploited to determine optimality. The optimality of the combined problem could then be 
declared if each individual problem was optimal on its own and each solution state and 
control history matched that of the combined solution. This verification was not 
performed in this work due to lack of time. 
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B. MAXIMUM SCIENTIFIC MASS 


The maximum scientific mass cost function added further difficulties in obtaining 
a solution. Simply changing the cost function yielded infeasible solutions. Unlike the 
minimum total fuel case, increasing the vehicle thrust did not resolve this issue. To 
obtain a solution, a hypothetical target planet with the properties of Mars was placed at 
5.2 AU. Like the minimum total fuel case, the problem was solved for a non-rotating 
target atmosphere. Even with these changes, solutions for high numbers of nodes became 
numerically unstable and yielded infeasible solutions. The solution presented below was 
obtained using 120 nodes for the low thrust segment and 50 nodes for the aerocapture 
segment. The state histories for this solution are given in Figure 65. The vehicle begins 
with a zero C 3 and initial mass of 660 kg and arrives at Mars with 487 kg consuming 173 
kg of propellant. The total trip time is 1187.4 days. The DIDO low thrust state histories 
are in excellent agreement with the propagated states. 


Low Thrust States: Earth to 5 AU 



Figure 65: Low Thrust State Histories (Maximum Scientific Mass) 
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Figure 66 shows the heliocentric trajectory for the low thrust portion. Again, the 
vehicle conducts a long continuous burn for 175 days before shutting off and coasting 
with just enough energy to reach 5.2 AU. 


Low Thrust: Earth to 5 AU 



AU 

Figure 66: Heliocentric Low Thrust Transfer Orbit (Max Scientific Mass) 


Figure 67 shows the control histories for the low thrust segment of the mission. 
The lower plot shows the normalized thrust history with distinct thrust switch at 
approximately 175 days. The thrust angle is approximately zero during the thrusting 
portion of the trajectory. The remainder of the thrust angle history may be disregarded as 
it is of course meaningless when thrust magnitude is zero. 
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Figure 67: Low Thrust Control Histories (Max Scientific Mass) 


The aerocapture state histories are show in Figure 68 and Figure 69. The 
trajectory begins with an inertial Y-infinity at arrival of 5211 m/s and a flight path angle 
of arrival of -10.1 degrees. The fact that the initial Yinfinity is not zero demonstrates 
that the global solution is minimizing low thrust propellant at the expense of more 
efficient thermal energy dissipation during the aerocapture segment. The initial heading 
of -95.6 degrees is due to the non-rotating atmosphere which causes the cost function to 
be invariant with latitude and heading angle. The inertial delta-V for the pass is 3668 m/s 
and terminates in an orbit with an eccentricity of 0.04. Capture occurs at 1.87 minutes. 
Unlike the low thrust trajectories, the DIDO solutions do not correspond as well with the 
propagated states, particularly with the radius state. It is likely that increasing the number 
of nodes for this segment would lead to better convergence between the DIDO and 
propagated solutions. 
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Figure 68: Aerocapture Position States (Max Scientific Mass) 


Velocity History 
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Figure 69: Aerocapture Velocity States (Max Scientific Mass) 
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The bank angle control history is shown in Figure 70 below. As in the pure 
aerocapture optimizations investigated earlier, the bank angle begins approximately lift 
up before switching at approximately 1.9 minutes to lift down. Again, the chatter in the 
early history of the control is due to the lack of aerodynamic control authority high in the 
atmosphere and would be expected to be smoothed with higher-node solutions. 


Control History 



Figure 70: Bank Angle History (Max Scientific Mass) 


Figure 71 shows the heating rate history during the atmospheric pass. The 
maximum stagnation point heating of 34.2 W/cm 2 occurs at 1.24 minutes and the total 
heat load is 3586.59 J/cm 2 . Using the TPS mass model discussed previously, this 
corresponds to a front-shield mass of 17.93 kg. Combined with the 9.83 kg of propellant 
required to circularize the orbit, the total aerocapture mass is 27.7 kg. This compares 
with 331 kg that would be required for a pure propulsive capture at the target planet. 
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Figure 71: Aerocapture Heating Rate (Max Scientific Mass) 
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XL CONCLUSIONS 


The suitability of the direct pseudospectral method for solving highly non-linear 
astrodynamic problems has been explored. For low thrust problems, the method has been 
shown to produce excellent results, particularly for minimum time problems. However 
unknown factors cause low thrust minimum fuel solutions to be more difficult to 
consistently obtain. The method was particularly successful in solving optimal 
aerocapture trajectories over a range of conditions. 

The suitability of the direct method for simultaneously solving a combined low 
thrust trajectory with terminal aerocapture was also explored. Although the fidelity of the 
models was reduced to obtain feasible solutions, the proof-of-concept is considered a 
success as it successfully found feasible solutions for the combined trajectories. This 
concept of simultaneously optimizing trajectory segments with vastly different dynamics 
has the potential to identify previously unexplored trajectory combinations and further 
research in this area is strongly suggested. 
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XH. FUTURE WORK 


Due to the proof-of-concept nature of this work, there are numerous areas for 
which future work is encouraged. For low thrust problems, launch vehicle optimization 
can be added allowing for design trades between initial mass and initial C3. More 
realistic trajectories may be obtained by exploring 3-DOF as well as non-circular target 
orbits. Additionally, adding the gravitational effects of additional bodies may allow for 
the exploitation of gravity assists and further mass savings. Finally, the difficulties with 
obtaining certain low thrust minimum fuel trajectories should be further investigated. 

For aerocapture, a more accurate TPS mass model should be developed such that 
both the heat load and the peak heating rate are taken into account. Furthermore, the 
benefits of angle of attack modulation during aerocapture as well as thrusting arcs may 
yield new families of trajectories and should be explored. Additionally, other cost 
functions such as minimum altitude may be useful for such missions as sample return. 

Finally, this initial work solving mixed dynamic optimization problems may be 
expanded in many areas. To obtain any solutions at all, the fidelity was considerably 
reduced. Hopefully many of the problems with the combined trajectories will be rectified 
when the difficulties with the minimum fuel trajectories discussed above are resolved. 
Lastly, with future versions of DIDO, users should be able to recover the covectors for 
problems with interior knots, allowing for the verification of the DIDO trajectories. 
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APPENDIX A: USEFUL TRANSFORMATIONS 


A. COORDINATE TRANSFORMATIONS 

The following transformations are useful for moving between various aerocapture 

frames. 

1. Spherical to Inertial 

First rotate about K by 0 , then Y' by -<|). 


e r 


COS(|) 

0 

— sin (f) 

COS0 

sin0 

O' 



= 

0 

1 

0 

-sin0 

COS0 

0 

J 

e. 

L <t> J 


sintf) 

0 

COS(|) 

0 

0 

1 

K 


which simplifies to 


COS0 COS(|) 

sin0 cost]) 

sincf) 


-sin0 

COS0 

0 

J 

-cos0 sincf) 

-sin0 sin ([) 

cost]) 

K 


2. REN Frame to Frenet Frame 

First rotate about r by\|/ , then rotate about cf)' by —y . 


e n 


cosy 

-siny 

O' 

'1 

0 

0 



= 

siny 

cosy 

0 

0 

cost)/ 

sint)/ 


_V 


0 

0 

1 

0 

-sint)/ 

cost)/ 

e, 

L <> J 


which simp lifies to 

cosy 
siny 
0 


-siny cost)/ 

-siny sin\)/ 


cosy cost)/ 

cosy sin \\r 


-sin \\i 

cost)/ 

e. 

L <I>J 



which can be inverted to 



cosy 

siny 

0 


-siny cost)/ 

cosy cost)/ 

-sint)/ 


-siny sint)/ 

cosy sint)/ 

cost)/ 

_K_ 


(213) 


(214) 


(215) 


(216) 


( 217 ) 
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which can be inverted to 



COS0 cost]) 
sin0 cost]) 
sin tf) 


-sin0 

-cos0 sintf) 


COS0 

-sin0 sintf) 


0 

cost]) 

e. 

L <l> J 


( 218 ) 


3. PCF to PCI Velocity 

Recall that Ypa is the inertial velocity vector and Ypcf is the velocity vector in 
the rotating frame. The transport theorem states 


V —V + PCI 9 PCF xr 
v pci v pcf + • xr 


But the angular rate of the rotating frame is given by 


PCI ? PLF = Q.K = sin tf) e r + cost]) & 


The positions vector in spherical coordinates is simply 

{r] = r <?„ 

The can be transformed from REN to spherical by 

{VpcU = 


cosy 

siny 

0 

" 0 " 

- sin y cost)/ 

cosy cost)/ 

-sin\|/ 

V 

-siny sin\|/ 

cosy sint|/ 

cost]/ 

_ 0 _ 


which gives 


{V PCF } = vsiny e r + vcosy cost)/ e e + vcosysint|/ e ^ 


so substituting into Eqn. (219) results in 

{Vpc/lrftj, 


vsiny 


e 

r 

4 

e, 

<i> 

vcosy cost)/ 

+ 

Q sintf) 

0 

£2cOSt|) 

vcosy sin\|/ 


r 

0 

0 


(219) 


( 220 ) 


( 221 ) 


( 222 ) 


(223) 


(224) 


which simplifies to 

{v,L = vsiny e r + (vcosy cost)/ +/T2 cost])) e 6 +v cosy sin\\r e ^ (225) 
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This gives the spherical inertial velocity components V r , V e , and V 0 in terms of 
the velocity vector in the rotating REN frame. 

R} ra =vsiny 

{Vg } ra =v cosy cost)/ +r£lc os(f) (226) 

K} ra =vcosy sinit)/ 


We can now obtain the velocity components in the inertial REN frame by squaring both 
sides of Eqn. (225) to get 

||y pc/ |r = v 2 + r'Q. 2 cos 2 cj) + 2rv£2cos(J) cost)/ cosy (227) 

so the inertial speed in the REN frame is 

V PCI = yjv 2 + r 2 Cl 2 c os 2 (|) +2 rv£2 cost]) cost)/cosy (228) 


Recalling Figure 3 


'E = atan 




(229) 


and 


f 


r = atan 


V 




(230) 


Substituting Eqns. (226) into (229) and (230) we get 


( 


'E = atan 


vcosy sin\)/ 


vcosy cost)/ +r£2cos(J> 


(231) 


and 


r = atan 


vsiny 


Jcos 2 y (v 2 +fl 2 r 2 )+ 2vcQcosy cost)/ cos (|) 


(232) 
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Thus Eqns (228), (231), and (232) give us the inertial velocity components 
resolved in the REN frame as functions of the position and velocity components as 
measured in the rotating REN frame. 
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APPENDIX B: MINIMUM AEROCAPTURE MASS AT MARS 
FROM ZERO ARRIVAL V-INFINITY DATA 


Data Summary: 

Inertial Velocity Components: 

Arrival V infinity (m/s): 0.00 

Atmospheric Entry: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg) 

Atmospheric Exit: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg) 

Delta-V (m/s) : 

Rotating Velocity Components 
Atmospheric Entry: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg): 


Atmospheric Exit: 


Speed (m/s): 

3274.53 

Heading (deg): 

0.09 

Flight Path Angle (deg): 

2 . 02 

Delta-V (m/s) : 

1427.83 


4702.36 

-0.09 

-7.83 


3523.53 

0.08 

1.87 

1425.77 


4949.30 

-0.09 

-7.44 


Trajectory Analysis: 
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Total Pass Time: 


734.40 sec (12.24 min) 


Minimum Altitude (km): 

Time to Min Alt: 

Max Dynamic Pressure (Pa): 

Time of Max Dynamic Pressure: 
Max Acceleration: (g) 

Time of Max Acceleration: 

Max Heating Rate: (W/cm A 2) 

Time of Max Heating Rate: 

Heat Load: (J/cm A 2) 

Final Orbit Parameters: 
Semi-major Axis (km): 

Periapsis (km): 

Apoapsis (km): 

Apoapsis Altitude (km): 
Eccentricity: 

Capture Time: 

Total Aerocapture Mass (kg): 
Propellant Mass to Circularize 
Estimated Front-shield mass (ki 


70.14 

158.94 sec (2.65 min) 

505.71 

158.94 sec (2.65 min) 

0.85 

158.94 sec (2.65 min) 

7.44 

140.18 sec (2.34 min) 

1700.38 

3563.65 

3437.39 
3689.92 
300.00 
0.04 

0.00 sec (0.00 min) 

19.14 

After Aerocapture (kg): 10.63 

): 8.50 


Pure Propulsive insertion: 

Isp = 330.000000 

Propellant Mass (kg): 201.442991 
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APPENDIX C: MINIMUM AEROCAPTURE MASS AT MARS 
FROM EXCESS ARRIVAL V-INFINITY DATA 


Data Summary: 

Inertial Velocity Components: 

Arrival V infinity (m/s): 

Atmospheric Entry: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg): 

Atmospheric Exit: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg): 
Delta-V (m/s): 

Rotating Velocity Components: 
Atmospheric Entry: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg): 

Atmospheric Exit: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg): 
Delta-V (m/s) : 


5200.00 

7178.83 

-0.15 

-10.47 

3517.05 
0.27 
2.18 
3661.78 


6933.98 

-0.15 

-10.85 

3268.11 

0.29 

2.35 

3665.88 


Trajectory Analysis: 
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Total Pass Time: 

Minimum Altitude (km): 

Time to Min Alt: 

Max Dynamic Pressure (Pa): 

Time of Max Dynamic Pressure: 
Max Acceleration: (g) 

Time of Max Acceleration: 

Max Heating Rate: (W/cm A 2) 

Time of Max Heating Rate: 

Heat Load: (J/cm A 2) 

Final Orbit Parameters: 
Semi-major Axis (km): 

Periapsis (km): 

Apoapsis (km): 

Apoapsis Altitude (km) : 
Eccentricity: 

Capture Time: 

Total Aerocapture Mass (kg): 
Propellant Mass to Circularize 
Estimated Front-shield mass (k< 


601.02 sec (10.02 min) 

58.51 

94.63 sec (1.58 min) 

2703.70 

79.78 sec (1.33 min) 

4.57 

79.78 sec (1.33 min) 

35.02 

72.76 sec (1.21 min) 

3388.37 

3550.25 

3410.57 
3689.92 
300.00 
0.04 

110.50 sec (1.84 min) 

28.75 

After Aerocapture (kg): 11.81 

): 16.94 


Pure Propulsive insertion: 

Isp = 330.000000 

Propellant Mass (kg): 386.401888 
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APPENDIX D: MINIMUM AEROCAPTURE MASS AT MARS 
FROM EXCESS ARRIVAL V-INFINITY DATA 


Data Summary: 

Inertial Velocity Components: 

Arrival V infinity (m/s): 

Atmospheric Entry: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg): 

Atmospheric Exit: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg): 
Delta-V (m/s) : 

Rotating Velocity Components: 
Atmospheric Entry: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg): 

Atmospheric Exit: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg): 
Delta-V (m/s) : 


9419. 99 

24987.43 

-0.08 

-9.35 

16143.37 

0.54 

1.28 

8844.06 


22615.43 

-0.08 

-10.34 

13736.80 

0.64 

1.51 

8878.62 


Trajectory Analysis: 
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Total Pass Time: 

1358.47 sec 

(22. 

64 min) 

Minimum Altitude (km) : 

256.46 



Time to Min Alt: 

249.75 sec 

(4.16 

min) 

Max Dynamic Pressure (Pa): 

3451.29 



Time of Max Dynamic Pressure: 

231.55 sec 

(3.86 

min) 

Max Acceleration: (g) 

5.83 



Time of Max Acceleration: 

231.55 sec 

(3.86 

min) 

Max Heating Rate: (W/cm A 2) 

264.59 



Time of Max Heating Rate: 

213.90 sec 

(3.56 

min) 

Heat Load: (J/cm A 2) 

55903.98 



Final Orbit Parameters: 




Semi-major Axis (km): 

24757.55 



Periapsis (km): 

23891.11 



Apoapsis (km): 

25624.00 



Apoapsis Altitude (km): 

1000.00 



Eccentricity: 

0.03 



Capture Time: 

213.90 sec 

(3.56 

min) 


327.98 
48.40 
279.52 

Pure Propulsive insertion: 

Isp = 330.000000 

Propellant Mass (kg): 528.735021 


Total Aerocapture Mass (kg): 

Propellant Mass to Circularize After Aerocapture (kg): 
Estimated Front-shield mass (kg): 
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APPENDIX D: MINIMUM AEROCAPTURE MASS AT MARS 

WITH G-LIMITS 

Data Summary: 

Inertial Velocity Components: 

Arrival V infinity (m/s): 5199.99 

Atmospheric Entry: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg) 

Atmospheric Exit: 

Speed (m/s): 

Heading (deg) : 

Flight Path Angle (deg) 

Delta-V (m/s) : 

Rotating Velocity Components 
Atmospheric Entry: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg): 


Atmospheric Exit: 


Speed (m/s): 

3271.83 

Heading (deg): 

-0.32 

Flight Path Angle (deg): 

2.16 

Delta-V (m/s) : 

3661.56 

Trajectory Analysis: 



6933.39 

-0.09 

-10.07 


3520.80 

-0.29 

2.01 

3658.02 


7178.82 

-0.09 

-9.73 


Total Pass Time: 


684.58 sec (11.41 min) 
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Minimum Altitude (km): 

Time to Min Alt: 

Max Dynamic Pressure (Pa): 

Time of Max Dynamic Pressure: 
Max Acceleration: (g) 

Time of Max Acceleration: 

Max Heating Rate: (W/cm A 2) 

Time of Max Heating Rate: 

Heat Load: (J/cm A 2) 

Final Orbit Parameters: 
Semi-major Axis (km): 

Periapsis (km): 

Apoapsis (km): 

Apoapsis Altitude (km) : 
Eccentricity: 

Capture Time: 

Total Aerocapture Mass (kg): 
Propellant Mass to Circularize 
Estimated Front-shield mass (k< 


61.98 

122.65 sec (2.04 min) 

1803.97 

96.30 sec (1.61 min) 

3.05 

96.30 sec (1.61 min) 

30.39 

78.32 sec (1.31 min) 

3814.61 

3557.99 

3426.09 

3689.90 

299.98 
0.04 

136.75 sec (2.28 min) 

30.20 

After Aerocapture (kg): 11.13 

r) : 19.07 


Pure Propulsive insertion: 

Isp = 330.000000 

Propellant Mass (kg): 386.401510 


144 



APPENDIX E: MAXIMUM INITIAL FLIGHT PATH ANGLE FOR 

AEROCAPTURE AT MARS 

Data Summary: 

Inertial Velocity Components: 

Arrival V infinity (m/s): 5200.00 

Atmospheric Entry: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg) 

Atmospheric Exit: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg) 

Delta-V (m/s) : 

Rotating Velocity Components 
Atmospheric Entry: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg): 

Atmospheric Exit: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg) 

Delta-V (m/s) : 

Trajectory Analysis: 

Total Pass Time: 708.99 sec (11.82 min) 


3436.81 

-69.58 

2.26 

3678.72 


7115.54 

-75.32 

-8.73 


3516.46 

-66.33 

2.21 

3662.37 


7178.83 

-73.46 

-8.65 
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Minimum Altitude (km): 

Time to Min Alt: 

Max Dynamic Pressure (Pa): 
Time of Max Dynamic Pressure: 
Max Acceleration: (g) 

Time of Max Acceleration: 

Max Heating Rate: (W/cm A 2) 
Time of Max Heating Rate: 

Heat Load: (J/cm A 2) 

Final Orbit Parameters: 
Semi-major Axis (km): 
Periapsis (km): 

Apoapsis (km): 

Apoapsis Altitude (km) : 
Eccentricity: 

Capture Time: 

Final Inclination: 

Total Aerocapture Mass (kg): 
Propellant Mass to Circulariz 
Estimated Front-shield mass ( 


66.52 

149.76 sec (2.50 min) 

1354.90 

114.93 sec (1.92 min) 

2.29 

114.93 sec (1.92 min) 

27.67 

91.29 sec (1.52 min) 

4612.22 

3549.02 

3408.13 

3689.92 

300.00 

0.04 

178.05 sec (2.97 min) 

70.00 deg 

34.97 

After Aerocapture (kg): 11.91 

r) : 23.06 


Pure Propulsive insertion: 

Isp = 330.000000 

Propellant Mass (kg): 386.401888 
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APPENDIX F: MINIMUM INITIAL FLIGHT PATH ANGLE FOR 

AEROCAPTURE AT MARS 

Data Summary: 

Inertial Velocity Components: 

Arrival V infinity (m/s): 5200.00 

Atmospheric Entry: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg) 

Atmospheric Exit: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg) 

Delta-V (m/s): 

Rotating Velocity Components 
Atmospheric Entry: 

Speed (m/s): 

Heading (deg): 

Flight Path Angle (deg): 


Atmospheric Exit: 


Speed (m/s): 

3358.39 

Heading (deg): 

-73.83 

Flight Path Angle (deg): 

4.71 

Delta-V (m/s) : 

3757.48 

Trajectory Analysis: 



7115.88 

-75.27 

-10.93 


3435.29 

-69.85 

4.60 

3743.54 


7178.83 

-73.40 

-10.83 


Total Pass Time: 


352.25 sec (5.87 min) 
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Minimum Altitude (km): 

Time to Min Alt: 

Max Dynamic Pressure (Pa): 

Time of Max Dynamic Pressure: 
Max Acceleration: (g) 

Time of Max Acceleration: 

Max Heating Rate: (W/cm A 2) 

Time of Max Heating Rate: 

Heat Load: (J/cm A 2) 

Final Orbit Parameters: 
Semi-major Axis (km): 

Periapsis (km): 

Apoapsis (km): 

Apoapsis Altitude (km): 
Eccentricity: 

Capture Time: 

Final Inclination: 

Total Aerocapture Mass (kg): 
Propellant Mass to Circularize 
Estimated Front-shield mass (k 


56.32 

93.33 sec (1.56 min) 

3298.44 

78.99 sec (1.32 min) 

5.57 

78.99 sec (1.32 min) 

39.76 

65.53 sec (1.09 min) 

3138.83 

3391.27 
3092.61 
3689.92 
300.00 
0.09 

98.27 sec (1.64 min) 

70.00 deg 

42.1 

After Aerocapture (kg): 26.41 

): 15.69 


Pure Propulsive insertion: 

Isp = 330.000000 

Propellant Mass (kg): 386.401889 
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